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A FINITE-DIFFERENCE PROGRAM FOR STRESSES IN 
ANISOTROPIC, LAYERED PLATES IN BENDING 

INTRODUCTION 


One critical feature associated with structural composites of laminated 
construction, using materials or geometrical arrangements that exhibit different elastic 
properties from layer to layer, is the possibility that the glued layers will separate or 
delaminate. This was undoubtedly realized from the outset of their use, and a brief 
historical sketch of the American scene is presented by Pipes [ 1 ] . However, the earliest 
serious investigation into the cause of delamination-type failure, namely the interlaminar 
stress problem, was apparently done in Japan by Hayashi [2,3], who reported that the 
maximum interlaminar shearing stresses occurred at the free edge of a laminate under 
tension. Hayashi used a plane stress model for the layers and approximated the 
interlaminar shears by a strain-averaging technique. Using a similar model, Puppo and 
Evensen [4] likewise discovered a sharp rise in the interlaminar stresses near a free edge. 
Notably, the use of the above models ignored the interlaminar normal stress. In two 
publications, Pipes and Pagano [5,6] developed a finite-difference program to solve the 
exact elasticity equations for a long laminate in uniaxial extension. In their development, 
the stresses are assumed independent of the axial coordinate and include all six 
components. The results of this investigation show that a sharp rise in both the 
interlaminar shear stresses and the normal stress occurs near the free edge. Thereafter, 
Oplinger [7] did an analysis of angle ply laminates in tension using a model similar to 
that of References 2 through 4. His approach allows a large number of layers to be 
considered. Indeed it was discovered that a singularity in the interlaminar shear occurs at 
the free edge of a laminate of one particular type of construction. An alternative solution 
to that employed in the above references is used by Rybicki [8] who applied a 
three-dimensional finite element formulation. His results agree with References 5 and 6. 

The present report marks the initial phase of a study of the interlaminar stresses 
induced in a layered laminate by bending. Following the approach used by Pipes [5], the 
laminate is modeled as a continuum and the resulting elasticity equations are solved using 
the finite-difference method. This solution technique is made possible by assuming that 
the laminate is bent into a cylindrical surface such that the stresses are independent of 
the axial coordinate. The objective of this report is to set forth the mathematical 
framework, present some preliminary results, and to avail the computer program to 
others. The results reveal a simplifying symmetry relationship in the displacements that 
will allow significant reduction in the size of certain numerical problems. In addition, 
trends in the interlaminar stress distribution are somewhat similar to those found for 
stretching problems, in that a sharp rise occurs at the free edge. 



PROBLEM FORMULATION 
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Laminate Description 

The laminated composites considered in this report consist of rectangular laminae 
symmetrically stacked with respect to a midplane and bonded together to form a flat 
laminate. The bonding is assumed to provide perfect adhesion between the laminae, 
which nullifies the possibility of slip between adjacent laminae thus establishing the 
conditions of continuous displacements and tractions at each interface. Each individual 
lamina is considered to be elastic, homogeneous, and orthotropic (i.e., each lamina 
possesses three planes of elastic symmetry). The assumption of homogeneity eliminates 
micromechanical effects such as those involving fibers or matrix. The geometry of a 
typical lamina and laminate is illustrated in Figure 1. One may note that the orthotropic 
coordinate axes (x',y',z) of a lamina are referred through a clockwise rotation about z to 
the fixed coordinate axes (x, y, z) of the laminate. The laminae are stacked along z to 
form a laminate whose sides are normal to x, y, and z. Each lamina is given a layer 
number m. 

/ 

Limiting the analysis to linear elastic materials, the constitutive relation for each 
lamina referred to the x, y, z coordinate system is 
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( 1 ) 


where the elastic constants Cjj are related to the nine orthotropic constants cjj through 

the well known transformation equations of References 9 and 10. 1 By associating the 
displacements u, v, and w with x, y, and z, respectively, the strains for each lamina are 
defined as 


1. In using the transformation equations in References 9 and 10 substitute -d for +6 since 
here the constants are referred to the unprimed coordinate axes of the laminate. 
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, ( 2 ) 


p m _ ,, m 

X ’X 


p m _ v m 

c y v >y 


e z m = w, z m 


Tyz m = w )y m + V, z m 7 xz m = w, x m + u, z m 7 


m = v m + m 
xy >x u >y 


where the comma denotes partial differentiation. 


Loading and Field Quantities 

Consider a laminate loaded by bending about y at the ends x = constant. 
Assuming that the laminate is long enough in the x-direction and that Saint-Venant’s 
principle holds for a laminate, the resulting stress distribution will be independent of x in 
regions sufficiently removed from the areas of loading. Using this assumption and 
following Lekhnitskii [11], the elastic strain-stress relations can be integrated to yield 
displacements for each lamina of the form 


u m = (C lY + C 2 z + C 3 )x + U m (y, z) 
v m = C,x 2 + C 4 xz + V m (y, z) 

w m = - X - C 2 x 2 - C 4 xy + W m (y, z) , (3) 

where U m , V m , and W m are unknown functions of y, z. The layer number, m, is left off 
the constants Cj because it results that each Cj must be the same for every lamina in 

order to satisfy the displacement continuity conditions at the interfaces. Thus, the Cj are 

found to be properties of the entire laminate. The displacement equations (3) represent 
the full three-dimensional elasticity solution that holds for all points in the laminate. 

To evaluate the Cj, the scheme is as follows. Since equations (3) hold for all 

points in the laminate, they must converge to the plane stress solution, which is an exact 
solution, in the interior region of the laminate. Integrating the relation [10,12] 


e i = B ii M j 


+ 


zD-jMj 


i, j = 1, 2, 6 


(4a) 
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for the case where Mj = -M and M 2 = M 6 = 0, the plane stress displacements are found to 
be 

Up S = (-Dj , Mz - B' 11 M)x - BgjMy - ^DUMyz + f(z) 
v ps = -^D' 16 Mxz - (B 21 M + D' 12 Mz)y + g(z) 

w ps = ^D' n Mx 2 + |D' 16 Mxy + ^D' 12 My 2 + f*(x) + g*(y) , (4b) 

where Bjj and D-j are laminate properties defined in Appendix A, and M is the applied 
moment. Comparing equations (3) and (4b) leads to the results: 

C x = 0 C 2 = -D'i i M 

(5) 

C 3 = -Bi,M C 4 = -|d' 16 M 

and 

U m (y, z) -» B u y + C 4 yz + U m (y, z) 

V m (y, z) -> B v y + D v yz + V m (y, z) 

W m (y, z) - -^D v y 2 + W m (y, z) , (6) 

where 2 

B u = -Bi.M , B y - -B 21 M , and D y = -D' 12 M . (7) 

2. The extended forms (6) for U m , V m , and W m are not necessary to the solution. 
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Substituting the results (6) into equations (3) yields displacements of the following 
functional form for each layer 


u m = (C 2 z + C 3 )x + (B u + C 4 z)y + U m (y, z) 
v m = C 4 xz + (B y + D y z)y + V m (y, z) 

w m = -|c 2 x 2 - C 4 xy - "D v y 2 + W m (y, z) , (8) 

where Cj, Bj, and D y are defined by equations (5) and (7). The strains are found by 

substituting the displacements (8) into the strain relations (2). The stresses then follow 
directly using the constitutive relation (1). 

It is of interest to examine the strain e x m which is 


e x m = C 2 z + C 3 . (9) 

Should the laminate be a balanced composite, i.e., the laminae are symmetrically stacked, 
according to composition and orientation with respect to the midplane z = 0, then B^j = 

0 and from equations (5) C 3 = 0, which results in a case of pure bending. For the 
opposite case, an unbalanced composite exhibits an extensional strain, C 3 , in bending. 
Such coupling effects are common to laminated composites. 


Field Equations and Boundary Conditions 

In regions sufficiently removed from the load planes, the nonboundary points 
must satisfy the reduced equilibrium equations 


m , m 

T i T 

xy,y 'xz,z 


a m + T m 

u y,y yz,z 


m 


r yz,y + °z,z 0 


( 10 ) 
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where the stresses exhibit no x-dependence, which conforms to an earlier assumption. 
Substituting for the stresses in terms of displacements yields the field equations for each 
lamina 


c™U™ y + cf 5 U™ + c™V™ y + c™V™ + (c™ + c™)W™ z = 0 


c™U™ y + c™U™ z + c™V™ y + cf 4 V™ z + (c™ + c™)W™ z = 0 


r m 4. „mu T m , /^m + .m w m , m w m , m w m 

(C36 c 4 s)Ujyz (C23 + C 44l’', yz T C 4 4 W 5 yy c 33 ”izz 


= -(c™C 2 + cf 3 D v + 2c™C 4 ) 


(ID 


The boundary conditions on the free surfaces normal to y are 


m - _m = m 


in = in _ 0 
' xy ' yz 


( 12 ) 


and on the free surfaces normal to z are 


a m = r m = r m = 0 
u z xz 'yz 


(13) 


For continuity at the interfaces, the boundary conditions are: 


(u m , v m , w m ) = (u m+1 , v m+1 , w m+1 ) 


and 


(a 7 m , 


r m 

XZ’ 



(14) 


respectively. 

It is noted that the corner conditions are ambiguous in that there are five possible 
conditions out of which only three can be employed at any one time. The remaining two 
may or may not be satisfied by the solution. Thus, combinations may be tried until some 
satisfying results are achieved. 
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FINITE-DIFFERENCE SIMULATION 


Function Representation 

The mathematical basis for the finite-difference method is Taylor’s Series. 
Referring to Figure 2, the Taylor Series expansion for a function f at some point y, z 
about the point (or node) I, J is 


f(y, z) = f(I, J) + yf, y (I, J) + zf, z (I, J) 

+ “ y 2f ,yy(I, J) + |z 2 f, zz (I, J) + yzf, yz (I, J) + ... • (15) 

Thus, for the specific node 1-1, J, the expansion is 

f(l-l,j) = f(l, J) - hjf.y + U^Tyy - ... . (16) 

Writing similar expansions for the remaining seven points neighboring the node I, J and 
simultaneously solving expansions for the first and second derivatives yields the 
finite-difference approximations for these derivatives. All but the last of these 
expressions, given below, are taken from Forsythe and Wasow [13]. They are 

f (i, j) = j)i + fd, j) + o(h z ) 

y h ! + h 2 Lh 2 hi J hi h 2 

f, z (I, J) = ~ [f(I, J + 1) - f(I, J - 1) ] + 0(h 2 ) 

f (I, J) = l— f-f f(I + 1, J) + -f- f(I - 1, J)1 - u - 2 — f(I, J) + 0(h 2 ) 

yy h! + h 2 |_h 2 hj J hj h 2 

f , zz a j) = [fa, j + 1) + fa, j - 1) - 2f(i, j)] + o(h 2 ) 

f, a, j) = 1 — r fd + 1, j + 1) - fa - 1, j + 1) - fa + 1, j - 1) 

y 2h 3 (h! +h 2 ) L 

+ f(I - 1, J - 1)] + 0(h 2 ) , (17) 
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where h is an order of magnitude equal to h 1; h 2 , or h 3 . The difference equations (17) 
are “central” differences. 

At boundaries and interfaces it is convenient to use “forward” and “backward” 
differences. Such difference equations are one-sided in that they express a boundary 
point in terms of neighboring points interior to the boundary. For the present problem, 
only first derivatives are of concern. 

To derive such difference equations, expand two points, both lying on one side of 
the reference point I, J, by using equation (15) in conjunction with Figure 2. For 
example, a forward expansion yields 

f(I + l,J) = f(I, J) + h 2 f, y (I, J) + ^ h 2 2 f, yy (I, J) + 0(h 2 3 ) 

f(I + 2, J) = f(I, J) + 2h 2 f, y (I, J) + ^ (4h 2 2 ) f, yy (I, J) + 0(h 2 3 ) . (18) 

Subtracting one expression from the other to eliminate the second derivative leads to the 
difference equation for the first derivative. Thus, the forward differences are 

f, y (I, J) = ~ [4f(I + l,J) - 3f(I, J) - f(I + 2, J)] - 0(h 2 2 ) 

f, z (I, J) = -j~ [4f(I, J + 1) - 3f(I, J) - f(I, J + 2)] - 0(h 3 2 ) . (19) 

Similarly, the backward differences are 

f,y(I, J) = [3f(I, J) + f(I - 2, J) - 4f(I-l,J)]+ 0(hj 2 ) 

f, z (I, J) = [3f(I, J) + f(I, J - 2) - 4f(I, J - 1)] + 0(h 3 2 ) . (20) 

It should be pointed out that more simplified, but less accurate, forward and backward 
expressions can be written; however, the present application requires all the accuracy that 
it is possible to attain near the free boundaries. Thus, the higher order difference was 
chosen. In addition, this choice yields a magnitude of error equal to that found in 
equations (17). 
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Using the representations just obtained, equations (11) through (14) can be 
transformed into difference equations characterizing the problem. For example, the last 
equation in (11) becomes 


2h 3 (hi + h 2 ) 


1 (c™ + cR) [U(I + 1, J + 1) - U(I - 1, J + 1) - U(I + 1, J - 1) 


+ U(I - 1, J - 1)] + (c^3 + cf 4 ) [V(I + 1, J + 1) 


- V(I - 1, J + 1) - V(I + 1, J - 1) + V(I - 1, J - 1)] 


2hj m 


hi + h 2 


cf 4 [w(I + 1, J) + ^ W(I - 1, J)] 


+ ^ c^ 3 [W(I, J + 1) + W(I, J - 1)] 

h 


,m , hj h 2 m 


- 2(c 44 + “J-^cy^Wd, J) = -h,h 2 [cf 3 C 2 + <#‘3 D 
h 3 2 


jn 


+ 2c§“ 6 C 4 ] , (21) 

where the layer number, m, is left off U, V, and W since their location is determined by 
the node I, J. 


Developing the Matrix Equation 

In this section, the difference equations, like (21), are transformed into a linear 
matrix equation of the form 


[A] [X] = [B] 


( 22 ) 


where A is an N X N coefficient matrix (N being the number of unknowns or equations), 
X is the solution vector, and B is the load or input vector. To accomplish this, the three 
unknowns (U, V, and W) must be uniquely collapsed into the single unknown X so that 
at each node three unique equations in X will be created. For instance, let 
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U-*X(1) 

V -> X(2) • 

W X(3) 


at Node 1 


U -»• X(4) 
V -»• X(5) 
W -»■ X(6) 


at Node 2 . (23) 


It remains to generalize such a transformation for all nodes. 

It is convenient to follow Pipes [1] and his notation is adopted. If LAT is the 
number of nodes in one column along the vertical axis (LAminate Thickness direction), 
then the nodes, unknowns, and equations can be identified by a unique number in terms 
of the nodal position (I, J). If 


JJ1 = 3[LAT(I - 1) + J] - 2 


(24) 


then 


NODE - LAT( I - 1) + J 
U(I, J) = X(JJ1) 

V(I, J) = X(JJ1 + 1) 

W(I, J) = X(JJ1 + 2) 


and 


(25) 


Number the 1st equation: JJ1 
Number the 2nd equation: JJ1 + 1 

Number the 3rd equation: JJ1 + 2 . (26) 

Letting 1=1 and J = 1 , 2 consecutively generates the results in (23). 

Since the finite-difference equations involve unknowns at nodes neighboring the 
JJ1 node, it is necessary to develop transformation relations like (24) in order to number 
unknowns at these nodes as well. For example, using I, J as the reference node, a 
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transformation relation for an unknown at the node I - 1 , J + 1 is found by letting I -> I 
- 1 and J -*■ J + 1 in (24) and giving the result a unique name, for example JJ7. Thus, 

JJ7 = 3 [ LAT(I - 2) + J] + 1 . (27) 

Using Table 1, which identifies all the unknowns at nodes neighboring I, J, and following 
the above procedure yields the transformation relations that uniquely number each 
unknown. In summary, all of these transformations are 

JJ1 = 3*(LAT*I1 +J) - 2 
JJ2 = 3*(LAT*I2 + J) - 2 
JJ3 = 3*(LAT*I2 + J) - 5 
JJ4 = 3*(LAT*I + J) - 2 
JJ5 = 3*(LAT*I + J) + 1 
JJ6 = 3*(LAT*I1 + J) + 1 
JJ7 = 3*(LAT*I2 + J) + 1 
JJ8 = 3*(LAT*I1 + J) - 5 
JJ9 = 3*(LAT*I + J) - 5 
JJ10 = 3*(LAT*I1 + J) - 8 
JJ11 = 3*(LAT*(I + 1) + J) - 2 
JJ12 = 3*(LAT*I1 + J) + 4 

JJ13 = 3*(LAT*(I - 3) + J) - 2 , (28) 


where 


11 = I - 1 

12 = I - 2 


(29) 
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TABLE 1 . NODE IDENTIFICATION 


Node 

U 

V 

W 

I, J 

X(JJ1) 

X(JJ1 + 1) 

X(JJ1 + 2) 

I - 1, J 

X(JJ2) 

X(JJ2 + 1) 

X(JJ2 + 2) 

I - 1, J - 1 

X(JJ3) 

X(JJ3 + 1) 

X(JJ3 + 2) 

I + 1, J 

X(JJ4) 

X(JJ4 + 1) 

X(JJ4 + 2) 

I + 1, J + 1 

X(JJ5) 

X(JJ5 + 1) 

X(JJ5 + 2) 

I, J + 1 

X(JJ6) 

X(JJ6 + 1) 

X(JJ6 + 2) 

I - 1, J + 1 

X(JJ7) 

X(JJ7 + 1) 

X(JJ7 + 2) 

I, J - 1 

X(JJ8) 

X(JJ8 + 1) 

X(JJ8 + 2) 

I + 1, J - 1 

X(JJ9) 

X(JJ9 + 1) 

X(JJ9 + 2) 

I, J - 2 

X(JJ10) 

X(JJ10 + 1) 

X(JJ10 + 2) 

I + 2, J 

X(JJ11) 

X(JJ11 + 1) 

X(JJ11 + 2) 

I, J + 2 

X(JJ1 2) 

X(JJ1 2 + 1) 

X(JJ12 + 2) 

I - 2, J 

X(JJ 1 3) 

X(JJ1 3 + 1) 

X(JJ 1 3 + 2) 


Generation of the matrix equation (22) now remains. To do this, straightforward 
substitution for U, V, and W, using Table 1, into equations (11) through (14) yields the 
desired results in equation form. For example, equation (21) becomes 
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,u lh h^{ (c ™ + c ™) [X(JJ5) - X(JJ7) - X(JJ9) + X(JJ3)] 

2h3(h! + h 2 ) 1 

+ (cft 3 + eft) [X(JJ5 + 1) - X(JJ7 + 1) - X(JJ9 + 1) 

+ X(JJ3 + i)] } + , 21l i— eft [X(JJ4 + 2) + ^ X(JJ2 + 2)] 

1 hj + h 2 hj 

+ cft[X(JJ6 + 2) + X(JJ8 + 2)] 

h 3 2 

- 2(cft + ^ eft) X(JJ1 + 2) 
h 3 

= -hi h 2 [eft C 2 + eft D y + 2cft C 4 ] . (30) 

To assure non-zero diagonal terms in the A-matrix, an appropriate equation number for 
(30) is JQ2 (in this case there is only one possibility) where 

JQ2 = JJ1 + 2 . (31) 

Now, from equation (30), the only nonzero elements for the JQ2 row in the A-matrix 
are 

A(JQ2, JJ5) = A(JQ2, JJ3) = C 
A(JQ2, JJ7) = A(JQ2, JJ9) = -C 
A(JQ2, JJ5 + 1) = A(JQ2, JJ3 + 1) = D 

A(JQ2, JJ7 + 1) = A(JQ2, JJ9 + 1) = -D 

A(JQ2, JJ4 + 2) = 2h 1 cft/(h 1 + h 2 ) 

A(JQ2, JJ2 + 2) = (h 2 /hj ). ■ 2h 1 cft/(h 1 + h 2 ) 

A(JQ2, JJ6 + 2) = A(JQ2, JJ8 + 2) = hjhjcft/hj 2 

A(JQ2, JJ1 + 2) = -2(cft + h 1 h 2 cft/h 3 2 ) , (32) 
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where 


C — h x h 2 (c™ + C4^)/2h 3 (h 1 + h 2 ) 

D = hxhjCc™ + c^4)/2h3(h! + h 2 ) • (33) 

Note that the material constants c™ 4 and c™ 3 are non-zero ensuring a non-zero 
diagonal element A(JQ2, JJ1 + 2). In addition to this, the load vector is 

B(JQ2) = -hjh 2 [cf 3 C 2 +cf 3 D v + 2 c^ 6 C 4 ] . (34) 


Of course, these results only apply to node numbers where the third equilibrium 
equation in (11) holds. The computer program logically connects appropriate equations 
with each node. The matrix elements for the remaining equations (11) through (14) are 
generated in a similar fashion. 


The Mesh Simulation 

The continuum is to be simulated by a number of nodal points that form a 
finite-difference mesh. The mesh is distributed over a cross section of the laminate as 
shown in Figure 3. The mesh is defined by the following parameters: 

NLAY : the number of laminae 

LAT: the number of nodes along one column in the LAminate 

Thickness direction 

LAW: the number of nodes along one row in the LAminate Width 

direction 

FSW1: the first change in nodal spacing termed Fine Simulation 

Width 

K: magnification factor of the fine simulation width 

H: the fine simulation width 
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Given these parameters, the following parameters can be determined: 

INF(M): values of J at the upper INterFace of the mth layer 

FSW2: the second change in nodal spacing 

KH: the coarse simulation width (K = 1, 2, 3, ...) 

JQMAX = 3*LAT*LAW: the number of unknowns or equations 

IBW = 2*(3*LAT + 1): the half bandwidth 

NBAND = 2* IBW + 1 : the full band 

The bandwidth of the coefficient matrix is found by considering that the maximum 
number of nodes involved in the difference equations is three, as can be seen from 
expressions (19) and (20), and calculating the maximum number of consecutive elements 
on both sides of the diagonal to and including the last off-diagonal non-zero element. 

Selecting equations representing the conditions to be imposed at each node 
remains to be accomplished. Because of the arbitrariness of the corner conditions, a 
number of choices are possible. Those selected for this program are illustrated in Figure 
4. 


A user’s guide and a more detailed description of the computer program are 
presented in Appendix C. A program listing is provided also in Appendix C. 


RESULTS 


The results given below were obtained using a square mesh, magnification factor 
K = 1, of size (LAW, LAT) = (13, 9). A complete mesh description, taken from the 
program output, is displayed in Table 2. It is seen that these dimensions represent a beam 
rather than a plate. The program was run on an IBM 370 computer utilizing virtual 
storage. 


A single material having properties typical of a high modulus graphite-epoxy was 
chosen for the above mesh. Using standard notation, 


E 2 2 = 20.0 X 10 6 psi , r i2 = v 13 = r 23 = 0.21 

E 22 = E 33 = 2.1 X 10 6 psi 

Gj 2 = G 13 = G 23 = 0.85 X 10 6 psi 
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TABLE 2. MESH DESCRIPTION TAKEN FROM PROGRAM OUTPUT 


-***- UNiFGRM SENDING OF A LAMINATED PLATE - 

*** INPUT CUl *** 

__ Mj*TE« OF LAYERS IN CROSS SECT f f Nt~ Nt~A¥ _ ~ Ar~ 

MJVHR OF NUDES ON VERTICAL AX I Si LAT = t? 

ML A TER If NODES ON HORIZONTAL AXIS, L AO- •=- ■ P ■■ 

CHVVJE IN M F SH WIDTH (FSWl) ATI ■= 3 

CHA\GE IN MESH WIDTH (FSW25 ATI = 7 

MESH WIDTH MAGNIFICATION FACTOR, K = I 


— - - 

LAYFl NO. 

l 

IN TER FACE 

AT 

J = 

4 

— — 

t A yER HO. - 

? — 

---- - ENTER FACE- 

AT- 

J-- - 

7 


LAYER NO. 

3 

INTERFACE 

AT 

J = 

10 


LAYER no--. - - 

A; 

1 NT ERF ACHE- 

TCP 

d 

13 


FINE- SI *FIL ATI ON WIDTH, ft = D.0016T-- 

where the subscript “1” refers to the fiber direction. The two laminate configurations 
which are considered are 

A = [ 6 , 0, 0, 9] 
and 

B = [0, 9, 9, 0] 
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with d as in Figure 1 such that 0 degree < d < 90 degrees. Typical laminate data and 
load constants are displayed in Table 3. 3 Here the additional constant MT is the resulting 
moment required to produce a specified maximum strain which, for the present analysis, 
is e x = 1.0 X 10~ 3 inch/inch (see Appendix B). 

A sample of the results for the displacement functions U, V, and W is presented 
in Table 4. Examination of their variation with respect to z reveals the apparent 
symmetry relations, 


U, V antisymmetric in z 

W symmetric in z 


within an accuracy of two digits. 

Symmetries with respect to y are evident for the strains within three-digit 
accuracy. Samples of these results are plotted in Figures 5 and 6. Coupling these apparent 
symmetries with the strain relations (2) in an expanded form yields 


U, V antisymmetric in y 

W symmetric in y 


The displacement results verify this precisely for U (to four places), but show some 
deviation in V and W. 4 

To illustrate the effect of bending on the stress distribution, Figures 7 through 19 
are presented. Although convergence to the exact values has yet to be demonstrated, the 
results do have qualitative merit. The following cases result from a bending strain of e x = 
1.0 X 10“ 3 inch/inch prescribed at the bottom surface. 

Of principal interest are the interlaminar stresses illustrated in Figures 7 through 
12. We note that laminates composed of 30 degree or 45 degree layers produce the 
greatest stress rise in <j z at the free edge with a more pronounced effect occurring if the 

angle plies are on the outside, i.e., system A = [0, 0, 0, 0] . A similar effect is seen in the 
shear stress Ty Z , although the rise in stress is sharply blunted by the requirement of zero 

3. The thermal problem is neglected in this preliminary analysis even though expansion 
coefficients appear in the program. 

4. It is interesting to note that the y-symmetries for V and W are verified precisely using the 
coarser mesh (LAW, LAT) = (8, 9) which decreases the relative size of the bandwidth. 
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TABLE 3. TYPICAL LAMINATE DATA AND LOAD CONSTANTS TAKEN 

FROM PROGRAM OUTPUT 


3A46P1M. GATA *** 


4.AVFR 

611 

C2? 

F33 

64-2 

-FI 3 

F23 

- HU1? 

NO 13 

+1023 

l 

20.GCCr+C> 

2. 103E+06 

2. lG3 r +36 

0.350F+06 

3. P506+96 

0. 850C f 06 

0.21 

0.21 

0.21 

2 

2C.0CC H * 36 

2. 10 3f + 06 

2. 1C0E + 06 

0. 853E+C6 

0. 350E+06 

0. 350F+06 

0.21 

0.21 

0.21 

3 

20.CCC5+36 

2.1 90E+06 

2.1 OOF +06 

0 .350F + 06 

3.P50E+06 

0. 850P+06 

0. 21 

0.21 

0.21 

4 

20.0C3E+35 

7.100E+06 

2. 1C0F+36 

0.350E+C6 

3. 950F+06 

0. 8 5 O r + 06 

0.21 

0.21 

0.21 


*** STIFFNESS MATR T C c S *** 

-4-A Y C ' l /TUFT A X - V - l ** ATR - f - X X - Y -Z ~ P R IM E M &T -R- i - 


-6. 7456+34 5.3456+06 6,2 l OF +05 0.0 
7 A SC +06 c «210CiQ5 -Ot8— 


0.0 4.536F+G6 2. 0740+0 3 * * * 7 5.6480+05 5. 6480+05 0.0 0.0 0.0 

0.0 4.506T I 00 2t ~ 2~ 1 30 + 06 4.771P » P5 0.0 o^o 0 . r 


2.2136+06 0.0 0.0 4.3876+06 

8. 500F+05 0.0 0.0 


9.500005 1.0 


2.2130+06 0.0 0.0 - o.O- 

-8.5000+06 0.-0- - - 0-.-O-- 


--8 v 6000+0 5 0.0 


R -.- 8 O 0 G -+ 05 * - 


- Jv - O - 


2.7136+06 4. 77 IF +03 0.0 0.0 

2.2136+04 0.0 -0.0 


- 0 .- 0 — 

0.0 

0.0 


-2-. 0 2 4 0 +3 7 5 . 66 00 + 0 5 3 . 66 8 0 * 0 6 0.0 


- 9 - . -A- 


2. 2330+04 4.77 10+05 0.0 0.0- - 0.-0- 

2.2130+06-0.0 0.0 -O.O-- 

8.500 G +Q5 0.0 0. 0— 


3. 509 6+05 3.0 

9. 500P+95 


8.5000+05 "O.-Q- 

8 . 8000+05- 


3 2. 0266+37 5. >4-36+05 5.648F+05 0.9 9.0 3.0 

? . 2 i 3E +04 4.771F+05 0.3 3.3 3.0 

-?. 2 O 3 F -+-> 4 - 0 ^- 4 - - — 9 -rO W 6 

0.0 9.509F+05 3.0 3.0 


3.6006+06 9 * (- 


2.0760+37 5.4480+05 5.4480+05 0.0 0.0 - -0*0- 

2.2130+06 6.7710+05 0.0 0.0 0.0- 

7.21 30 + 06 0.0 OvO ■ A vO 

8.8000+06 9.0 0.0 

^ 000 + 0 ^ - 0 . 6 - 


). ^ 0 3 P + 35 


«. *00 r '+05 



TABLE 3. (Concluded) 



TABLE 4. DISPLACEMENT FUNCTION RESULTS TAKEN FROM 
PROGRAM OUTPUT FOR LAMINATE DESCRIBED IN 
TABLES 2 AND 3 


OESE^AEEi^ENE- FUNErTf OHS ***-- - - - 

\C 9E J-ni SPLACEMFNT V-PISPLACEMENT W-nrSPLACEMENT 


1 

0. 161 56 ID— 04 


0 .2646360- 34 

-0. 9090399-05 

7 

0. 1495800-04 


0.2198310-04 

-0.9368049-05 

3 

0. 1253810—04 


0.1731769-04 

-0.9519399-05 

\ 

0.9535940-05 


0 . 1270140-04 

-0. 9535900-05 

5 

0.61 16960-05 


0. 3189079- 35 

-0. 9400029-05 

A 

0. 3043950-05 


0.4033640-05 

-0. 9246869-05 

1 

0. 4871890-09 


0.2507690-03 

-0.9165899-05 

£ 

-0. 3042910-05 


-0 .4039139-05 

-0.925034^-05 

~3 

-0.6115750-05 


-0.3194320-05 

-0.9376939-05 

11 

-0. 9? 6689 0— 05 


-0 . 1 263080- 34 

-0. 9528929-05 

l 1 

-0. 1253540-04 


-0 .17321 19-04 

-0. 9544139-05 

t 2 

-0. 1495509-04 


-0.2193309-04 

-0.9371619-05 

l 3 

-0. 1615030-04 


-0.2648090-04 

-0.9093919-05 
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stress at the free edge, and here the stress in system B = [0, 0, 0, 0] is slightly more 
pronounced than that in A. The largest stress rise, an order of magnitude greater than o z 

and Ty Z , is created in the A-system in r xz . Again it is the 30 degree laminate incurring 

the sharpest stress rise, but here the 15 degree laminate overshadows the 45 degree 
laminate. In summary, the laminates containing 15 degree through 45 degree layers 
located adjacent to 0 degree layers have the largest interlaminar stresses for the cases 
considered; i.e., 0 degree < 0 < 90 degrees taken in 15 degree intervals. 

Some results peculiar to the numerical method of solution should be pointed out. 
Referring to Figure 9, we note a sharp rise in the stress ct z at the midpoint node (I, J) = 

(5, 7). This is a result of fixing the displacements at I = 5 and 6, J = 7 in the program in 
order to zero-out rigid body motion and drift in the solution routine. However averaging 
the values for ct z just above and just below the interface (at J = 7, m = 2 and m = 3) 

yields a more plausible result. Since the tractions must be continuous at the interface 
anyway, this averaging technique was also applied at the free edges where the free surface 
conditions were adopted in lieu of the continuity conditions. This technique had varying 
success as illustrated by the 75 degree and 90 degree configurations in Figures 10 and 11. 

The in-plane stresses are illustrated -in Figures 13 through 19. In Figure 13, we 
find that o x in the 0 degree layers is independent of the orientation of the adjacent layer 

when the maximum strain is specified. 5 This facilitates the presentation of both systems 
A and B in one figure. It is interesting to note in Figure 14 that a x rises at the free edge 

if the 0 degree layers are outside the laminate and drops if these layers are inside the 
laminate. 

Observation of Figures 15 and 17 for the distribution of a y and r xy with respect 

to z reveals that the off-axis layers, particularly again for 15 degrees through 45 degrees, 
serve as stress raisers with the effect considerably more pronounced if the 0 degree layers 
are inside. 

Typical distributions of a y and r xy with respect to y are shown in Figures 18 and 

19. The disturbing feature of these plots is that the stresses just above an interface do 
not approach zero at the free surface. One cause of this problem is the placement of 
nodes directly on the interface, which requires their occupation by both layers. Then at 
the corners, as stated previously, the multitude of boundary conditions cannot be 
satisfied. 6 However this problem is confined to the free surface nodes and one line of 


5. In agreement with the beam theory approximation. 

6. Placing the interface between two nodal lines may alleviate this problem. 
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interior nodes. To see this, one may examine the curves for the A-system at J = 4- and J 
= 10+ and note that they are reflections of each other within the range 3 < I < 7. Since, 
from above, Oy and r x y appear, in general, to be antisymmetric in z, the correct values at 

1=10+ are recovered within this range if we accept the values at J = 4-. 


CONCLUSIONS 


Although only two types of laminate systems were considered, namely A = [ 0 , 0, 
0, 0 ] and B = [0, 0 , 0 , 0], it is reasonable to assume from these results and from 
physical considerations that the following symmetry relations hold for balanced (Bjj = 0) 
composites: 


U, V antisymmetric in y and z 
W symmetric in y and z 


where U, V, and W are displacement functions of y and z. Based on the stress results, 
laminates containing layers oriented within the range 15 degrees < 0 < 45 degrees 
produce the largest interlaminar stresses out of the cases studied, 0 degrees < 0 < 90 
degrees taken in 15 degree intervals. In fact this same group of laminates produces high 
values in the in-plane stresses as well, with the effect considerably more pronounced for 
the A-system. Although some deviations in stress occur in the numerical solution, they 
are localized to a double line of nodes at the boundary. This is a disconcerting feature of 
the solution in that the boundary region stresses appear to be critically involved in 
delamination-type failure, which makes their accurate determination desirable. 

This study provides a base for future work in this area. Using the present program 
coupled with an out-of-core equation solver routine, unbalanced laminates may be 
studied. Using the symmetry relations discussed above, the present computer program 
may be modified to more efficiently handle balanced laminates (Bjj = 0). 
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Figure 2. 
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= 0 


*AT INF(m) WHERE 1 < I < LAW AND 1 < m < NLAY : 
[u m , v m . w m ] = [u m+1 , v m+1 ,w m+1 ] 


r m i 

[a . T . T 
1 z ' yz' 1 xz J 


[a m+1 r m+ 1 , T m+1 ] 
[a z ’ r yz ' xz 


I = l 


STATIC EQUILIBRIUM IS IMPOSED AT ALL INTERIOR POINTS 


Figure 4. Equations selected for each node. 
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Figure 8. Variation of the normal stress a z (symmetric in y) with y 
for a [0, 6, 6, 0] laminate. 
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Figure 10. Variation of the shear stress Ty Z (antisymmetric in y) 
with y for a [0, 0, 0, 0] laminate. 
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-1 0 1 2 3 ' ' 4 6 8 10 12 .14 16 

Oy (psiX 10 2 ) 

Figure 15. Variation of the normal stress o (antisymmetric in z) 
with z for a [0, 0, 0, 9] laminate. 












Figure 18. Variation of the normal stress Oy with y. 
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APPENDIX A 


LAMINATE CONSTANTS 


Following Reference 9 or 10, define 




c i 3 m c j3 m 


-33 


.m 


i, j = 1,2,6 


and let t be the half-thickness of the laminate, ho the thickness of a lamina, and N the 
total number of layers; then 


N 

nj = ll ° 22 QiJ 

m=l 


A;; = 


m 


h 2 

B-. = 

« 2 


f N 


N 


22 Qij m(2m - d - n 22 Qi 


m 


J. m=l 


m=l 


D- = 

y 3 


N 


22 Qij m (3m 2 - 3m + 1) 
m=l 


N 

- § N 22 Qij m (2m - 1) 

m=l 


♦ £ Q ij 

m=l 


m 


with i, j = 1, 2, 6. Finally, let 


A* = A' 1 , B* - -A" 1 B , and D* = D - BA" 1 B 
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where the letters symbolize 3X3 matrices. Then, 


B' = B*(D*) _1 
and 

D' = (D*)' 1 

Considerable simplification is attained if the laminate is balanced, which implies B 
= 0 . 



APPENDIX B 


STRAIN SPECIFICATION 


Rather than prescribe the laminate' loading as end moments, the maximum strain, 
e x max , at the top and bottom surfaces, z = ±z max , will be prescribed. From equation (9), 

we have 


max 


= Co z max + C 


Now from equations (5), 


C 3 = -BWM = ^C 2 

u i 1 


so that 


max _ 


c, i z max + 

u i i 


and, thus, 


C, = 


tv max 
U 11 e X 


B 1 1 + D 1 1 z 


' ^max 


In the computer program, we set e x max = -1.0 X 10' 3 inch/inch at the top surface z = 
+ z max to eva i ua te the constant C 2 which represents the inverse bending radius. 
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APPENDIX C 


THE COMPUTER PROGRAM 


Program Description 

The computer program is an in-core program and is not overlayed. It is felt that a 
flow chart of the program would be no less complicated than the presentation of a listing 
with an accompanying explanation, so the latter choice will be followed. Certain 
statements in the program are extraneous to the problem in this report because the 
program is in steady transition to handle more general problems. A part-by-part 
description follows. 

Part I. Part I contains a brief definition of terms and an explanation of the order 
and format of the data cards. The dimensions of the data are: H is in inches, E is 
material constants in psi (the shear moduli Gi 2 , etc., are read into the El 2, etc., arrays), 
ALPHA is the coefficient of expansion in inches/inch/°F, and THETA is the lamina 
orientation in angular degrees. Precision and dimension statements are then established, 
data are read in, and mesh parameters are' calculated. The letter M refers to the layer 
number. In the loop, DO 9000, IRAN counts each laminate layup from one to IRUN 
(only changes in lamina orientation are allowed for within this loop). 

Part II. Part II calculates the anisotropic stiffness matrix. BETA is in radians. 
CP11, etc., are the orthotropic elastic constants in the primed coordinate system. 
Cll(M), etc., are the anisotropic elastic constants for the Mth lamina in the x, y, z 
coordinate system. AL1P(M), etc., are the coefficients of thermal expansion in the 
primed coordinate system and AL1(M) are those coefficients in the x, y, z coordinate 
system, both the the Mth lamina. Finally, the subroutine MATCON, which calculates the 
laminate MATerial CONstants, is called. 

Part III. Part III calculates the coefficient matrix for the difference equations. 
The loops DO 100 and DO 101 count through the mesh node-by-node. DO 3000 zeroes 
out the A-matrix. 

The logic that associates the various field conditions with each node and correctly 
fills out the A-matrix is contained in DO 102. First the node I, J is tested to determine 
the proper layer number, M. Then the node is checked to see if it lies on a boundary, 
along J equals a constant, or lies at some select position (in this case, IMID or JMID). If 
it does, the program is routed to the statement number that contains the non-zero matrix 
elements satisfying the conditions imposed at this node. Should the node not lie at any 
of these preselected locations, the program passes through the IF statements on J to 
statement number 193, which initiates a series of checks to see if the node lines on 
selected values of I. These values include the boundaries I = 1 or I = LAW, the changes in 
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nodal spacing I = FSW1 or I = FSW2, and all points in the region between FSW1 and 
FSW2. Should the node not lie at any of these locations, the program passes through the 
IF statements on I and evaluates the non-zero coefficients for the only remaining 
possibility, the equilibrium terms for a square mesh. 

When a node does he on some select location, say J equals LAT, then the logic in 
that statement series, say the series starting from statement number 202, guides the 
program through the checks on selected values of I in a fashion similar to that above. 
The logic is easily understood by reading directly from the listing. 

Upon reaching statement number 102, the A-matrix (3 X JQMAX) is full. The 
elements of the A-matrix lying within the bandwidth are then stored in the banded 
matrix AX. The loops DO 100 or DO 101 then continue for the next node, if any. The 
previous A-matrix is destroyed and regenerated for the new node until the loop DO 100 
is satisfied. 

At rewind 9, the matrix AX and the load vector X are stored for later use. The 
loop DO 107 stores the load vector X(I) in AX(I, NBD). Then a series of WRITE 
statements (listed as comments) will output the coefficient matrix AX and load vector X 
should they be desired. Finally, the solver routine, TRMSTR, is called. 7 

Part IV. Part IV outputs the functional displacements and provides an accuracy 
check. Just below statement number 4006, the STOP 1 statement will terminate the 
program if the coefficient matrix AX is singular. (Such an occurrence probably indicates 
an error.) The loop DO 108 stores the solution vector AX(I, 1) in X(I). Then the original 
values of the matrix AX and load vector X are read back into the AX array and R 
vector, respectively. 

The loops DO 11 and DO 12 output the values for the functions U(y, z), V(y, z), 
and W(y, z) which occur in the displacements u, v, and w, respectively. 

The series of statements from the one above 9950 to 9990 outputs the accuracy 
results. These results provide the difference between the original load vector, now stored 
in the R-array, with the calculated load vector, which is found by substituting the 
appropriate solution vectors, X(I), into each matrix equation. In addition to giving the 
accuracy of each equation, an average accumulated accuracy is provided. 

Part V. Part V outputs the strains and stresses. The logic is similar to that in Part 
III. Knowing that the finite-difference relations for the strains differ for various mesh 
locations, the strains are split into terms dependent upon the value of I and terms 
dependent upon the value of J. The strain SX, which represents e. x , depends upon neither 

the value of I nor the value of J and is determined prior to any logical branching. 


7. Actually the AX-matrix stores a transposed A-matrix; i.e., instead of storing row elements 
crosswise or in a row, they are stored in the AX-matrix vertically or in a column. The result 
is a drastic reduction in “wall-time” on the IBM 370. This necessitated a slight revision in 
the solver routine, TRIMSS, as written by Billy Gibbs, U.S. Army, Redstone Arsenal [14]. 
So here it is called TRMSTR or TRIMSS transposed. 
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First, the node is checked to determine its location with respect to I, and 
I-dependent strains (or the partial strain, SYZI) are calculated. Then the loop DO 392 
establishes the correct layer number, M, in order to check if J lies on the interface, 
INF(M). Upon determining the correct location of the node with respect to J, the 
J-dependent strains (or the partial strain, SYZJ) are calculated. Statement number 391 
totals the partial strains to obtain SYZ. The stresses are then calculated in a straight 
forward manner using equation (1). Note that the stresses are calculated twice at 
interface nodes, once for the material below the interface and again for the material 
above. 


Part VI. The subroutine MATCON calculates the MATerial CONstants Cj, BU, 
BV, and DV as defined earlier in the text. 

Part VII. The subroutine MAMULT is a MAtrix MULTiplier and is easily 
understood from the listing. 

Part VIII. The subroutine MATIN4 is a MATrix INversion routine which is 
described in Reference 14. 

Part IX. The subroutine TRMSTR is the equation solver which is described in 
the listing. 

Part X. The subroutine RITE is used to wRITE out a matrix or vector. 


Program Listing 

The complete listing of the program is contained in the following pages. 
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FORTRAN IV G1 RELEASE 2.0 


MAIN 


DATE = 75007 


08/16/07 


0006 

0007 


0008 


0009 

0010 
0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 


00000010 
00000020 
00000030 
00000040 
00000050 

IF THE NUMBER OF LAYERS EXCEED 6, THE COMMON /MC / AND DIMENSION t E 1 1 *00000060 
E22* ETC.) STATEMENTS MUST BE REDIMENSIONED TO AGREE WITH LAT. 

REMEMBER TO PLACE A COMMON /MC/ STATEMENT IN SUBROUTINE MATCON. 


JQMAX IS THE NUMBER OF UNKNOWNS OR EQUATIONS TO BE SOLVED. 

A IS A FULL MATRIX (3 X JQMAX) REPRESENTING EACH NODE. 

AX IS THE BANDED MATRIX (NBAND+1 X JQMAX). 

X IS THE LOAD VECTOR. AFTER TRIMSS X BECOMES THE SOLUTION VECTOR. 


USE THE FOLLOWING ORDER FOR DATA CARDS 


DATA CARD 
1 
2 

3 

4 

5 

NOTE* 


NO. 


DATA 
LAT, LAW, FSW1, 


NL AY, 

H 

Ell* E22, E33, E12, E13* E23 
NU 12, NU 13, NU23 

ALPHA 1 PRIME, ALPHA 2 PRIME, ALPHA 3 PRIME 
REPEAT CARDS OF THE TYPE 3, 4, 5 FOR EACH ADDITIONAL 


00000070 
00000080 
00000090 
00000100 
00000110 
FORMAT 00000120 
5110 00000130 

G12.5 00000140 

8G12.5 00000150 
8G12.5 00000160 
8G12.5 00000170 
/ER 00000180 



C 

6 

SXMAX 

, C3E 

10G10.3 

00000190 


C 

7 

IRUN 


5110 

00000200 


C 

8 

THETA ( 1 ) , THETA ( 2 ) , THETA ( 3 ) , ETC. 

10G10.3 

00000210 


c 

NOTE, REPEAT CARD 8 

FOR EACH ADDITIONAL LAYUP. 


00000220 


c 





00000230 


c 





00000240 

0001 


INTEGER 

P, FSW1, 

FSW2 


00000250' 

0002 


DOUBLE 

PRECISION 

TEST, R, ERR, AVE, DT 


00000260 

0003 


DOUBLE 

PRECISION 

AX, X 


00000270 

0004 


DOUBLE 

PRECISION 

THETA, BETA 


00000280 

0005 


DOUBLE 

PRECISION 

CM, CN, CM4, CN4 , CM3N, CN3M , CM2, 

CN2 , GNU21, 

00000290 


GNU31 , 
CP23, 


GNU32 , DET, CPU, 
CP44, CP55, CP66 


CP22, CP33 , CP 12 , 


DIMENSION AX( 162,351) , A ( 3 , 35 1 ) » X{351), R(351> 

COMMON /MC/ Cll(6) , C12 < 6) , C16< 6 ) , C22 ( 6 ) ,C26 ( 6 > ,C66 ( 6 > , C 13 < 6 ) , 

1 C23I 6) , C36 ( 6 ) ,C44(6) ,C45(6) , C 55 ( 6 ) , C33 < 6 ) , AL 1 1 6 ) , AL2 I 6 ) , 

2 AL3(6) , AL6( 6) ,C2,C3, C3E , C4, BU , DU , BV , DV , H, SXMAX , NL AY , INF (6) 

DIMENSION E 1 1 < 6 ) ,E22(6),E33(6) ,E12(6),E13(6),E23(6), GNU 12 ( 6 ) , 

1 GNU13(6) ,GNU23{6) ,THETA(6) , AL1P16), AL2P(6), AL3P(6) 

TEMP = 0.0 

WR I TE ( 6 , 600 ) 

READ { 5 , 601 )NLAY,LAT»LAW,FSW1,K 

FSW2=LAW-FSW1+1 
JQMAX = 3*LAW*LAT 
IBW = 2*(3*LAT+1) 

IBW1 = IBW+1 
NBAND = 2*IBW+1 

WRITE! 6,602 )NL AY, LAT, LAW, FSW1,FSW2,K 

LAT1=LAT-1 

IMID = ( L AW + 1 ) / 2 

JMID = (LAT+1J/2 

DO 501 M= 1 , NL A Y 


CP13, 00000300 
00000310 
00000320 
00000330 
00000340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000400 
00000410 
00000420 
00000430 
00000440 
00000450 
00000460 
00000470 
00000480 
00000490 
00000500 
00000510 
00000520 
00000530 
00000540 
00000550 
00000560 
00000570 
00000580 
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0022 


INF ( M ) » 1+M*LAT 1/NLAY 


00000590 

0023 


WRITE(6,608)M, INF(M) 


00000600 

0024 


501 CONTINUE 


00000610 


C 



00000620 


c 

NOTE THAT INF(NLAY) EQUALS LAT AND IS NOT AN ACTUAL INTERFACE. 

00000630 


c 



00000640 

0025 


READ ( 5,603 ) H 


00000650 

0026 


WRITE ( 6, 607 1 H 


00000660 

0027 


HSQ = H**2 


00000670 


c 



00000680 

0028 


WR ITE ( 6 1 604 ) 


00000690 


c 



00000700 

0029 


DO 500 M*1 ,NLAY 


00000710 

0030 


READ(5,603)E11(M) ,E22(M) ,E33(M) ,E1Z<M) ,E13<M) ,E23(M) 


00000720 

0031 


READ! 5,603 )GNU12(M) , GNU 131 M) *GNU23(M) 


00000730 

0032 


WR I TE ( 6, 605 ) M, Ell(M), E22(M), E33(M), E12(M), E13(M), 

- E23 ( M ) t 

00000740 



1 GNU1 2 ( M ) , GNU 1 3 ( M ) , GNU23(M) 


00000750 

0033 


READ(5,603) AL1P{M) , AL2P(M>, AL3P(M) 


00000760 

0034 


500 CONTINUE 


00000770 


c 



00000780 

0035 


RE AD ( 5,606) SXMAX, C3E 


00000790 

0036 


R E AD ( 5 f 60 1 ) IRUN 


00000800 


c 



00000810 

0037 


DO 9000 IRAN = 1, IRUN 


00000820 

0038 


READ ( 5,606) ( THETA ( M ) ,M= 1 ,NL A Y ) 


00000830 


c 



00000840 


£***#* ***** ******** **************** ************ ** *****#*****##***#****#*000008 50 


C 



00000860 


C 

CALCULATION OF ANISOTROPIC STIFFNESS MATRIX TERMS REFERRED 

TO X » Y » Z 

00000870 


C 



00000880 


C- 

**######## #############*#*##*##*##*####*#*#**###*#####*## ###***##**####QQQQQg90 


C 



00000900 

0039 


WRITE(6,613) 


00000910 

0040 


XX = 0.0 


00000920 

0041 


DO 3001 M= 1 ,NLAY 


00000930 

0042 


BETA = .0174532925199433D0*THETA(M) 


00000940 

0043 


CM=DCOS ( BETA) 


00000950 

0044 


CN=DS I N ( BET A ) 


00000960 

0045 


IFIDABS(CM) . L T . 1 # E —0 8 ) CM = 0. 


00000970 

0046 


IF (DABS(CN) . LT. l.E-08) CN = 0. 


00000980 

0047 


CM4=CM**4 


00000990 

0048 


CN4=CN**4 


00001000 

0049 


CM3N=CM**3*CN 


00001010 

0050 


CN3M=CN**3*CM 


00001020 

0051 


CM2=CM**2 


00001030 

0052 


CN2=CN**2 


00001040 

0053 


GNU21=GNU12(M)*E22(M) /Ell (M) 


00001050 

0054 


GNU3 1 = GNU 1 3{M)*E33(M)/E11(M) 


00001060 

0055 


GNU32=GNU23(M)*E33(M)/E22(M) 


00001070 

0056 


DET = 1 . -GNU1 2 { M ) *GNU2 1— GNU23 (M)*GNU32-&NU13tM)*GNU31 


00001080 



1-2.*GNU12(M)*GNU23(M) *GNU31 


00001090 

0057 


CP11=E11(M)*( l.-GNU23(M) * GNU 3 2 ) /DE T 


00001100 

0058 


CP22=E22(M)*( l.-GNU13(M)*GNU31)/DET 


00001110 

0059 


C P33= E 33 ( M ) * ( l.-GNU12(M)*GNU21 J/DET 


0000 1120 

0060 


CP12=E11<M)*(GNU21+GNU23(M)*GNU31 )/DET 


00001130 

0061 


CP13 = E1UM)*(GNU31+GNU21*GNU32)/DET 


00001140 

0062 


CP23=E22(M)*(GNU32+GNU12(M)*GNU31 )/DET 


00001150 

0063 


C P44 = E 23 ( M ) 


00001160 
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0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 


0079 

0080 
0081 
0082 

0083 


0084 

0085 

0086 

0087 

0088 
0089 


0090 

0091 

0092 

0093 

0094 


IV G1 RELEASE 2.0 


MAIN 


DATE = 75007 


08/16/07 


CP55-E 13 { M ) 00001170 
CP66=E12 ( M ) 00001180 
C11(M)=CM4*CP11+2.*CM2*CN2*CP12+CN4*CP22+4.*CM2*CN2*CP66 00001190 
C12(M>=CM2*CN2*CP11+(CM4+CN4)*CP12+CM2*CN2*CP22-CM2*CN2*4.*CP66 000012.00 
C16<M)=CM3N*CP11~(CM3N-CN3M)*CP12-CN3M*CP22-2.*(CM3N-CN3M)*CP66 00001210 
C22(M)=CN4*CP11+2.*CM2*CN2*CP12+CM4*CP22+4.*CM2*CN2*CP66 00001220 
C26<M)=CN3M*CP11-(CN3M-CM3N)*CP12-CM3N*CP22-2.*(CN3M-CM3N)*CP66 00001230 
C66(M)=CM2*CN2*CP11-2.*CM2*CN2*CP12+CM?*CN2*CP22+(CM2-CN2)**2*CP6600001240 


C13(M)=CM2*CP13+CN2*CP23 
C231M)=CN2*CP13+CM2*CP23 
C36 ( M) =CM*CN*( CP 13-CP23 ) 
C44(M)=CM2*CP44+CN2*CP55 
C45 ( M)=CM*CN*( CP55-CP44) 
C55(M)=CN2*CP44+CM2*CP55 
C33 ( M ) =CP33 


00001250 

00001260 

00001270 

00001280 

00001290 

00001300 

00001310 

00001320 

00001330 


£**************************************************** ***************** **00001 340 
C 00001350 
C CALCULATION OF THE COEF. OF THERMAL EXPANSION REFERRED TO X,Y,Z 00001360 
C 00001370 

£***********************************************************************00001380 

C 00001390 

00001400 
00001410 
00001420 
00001430 
00001440 


AL1 (M)=CM2*AL1P(M}+CN2*AL2P(M) 
AL2(M)-CN2*AL1P(M)+CM2*AL2P(M) 
AL3(M)=AL3PIM) 

AL6(M)=2.*CM*CN*( AL1P< M)-AL2P( M) ) 


C36 { M > T CP33, XX, XX, XX, 
XX, CP44, XX, XX, C 55 { M ) , 


XX, 


C 16 ( 

M) , 

CPU, 

00001450 

M) , 

XX, 

XX, 

00001460 

M) , 

XX, 

XX, 

00001470 

C44( 

M) , 

C45 ( M ) , 

00001480 

XX, 

C66I 

[M), CP66 

00001490 


3001 CONTINUE 


00001500 
00001510 
00001520 
00001530 
00001540 
00001550 
00001560 
00001570 
00001580 
00001590 
00001600 
00001610 
00001620 

**********************************************************************00001630 

00001640 

CALCULATION OF THE COEFFICIENT MATRIX FOR THE DIFFERENCE EQUATIONS 00001650 

00001660 

**********************************************************************00001670 

00001680 


WR I TE { 6, 6 1 1 ) 

DO 503 M = 1 , NLAY 

WR I TE ( 6 , 614 ) M, THETA ( M ) , AL1(M), AL2 ( M ) » AL3 ( M ) , AL6 ( M ) , 
1 AL1P(M), AL2P(M), AL3P l M ) 

503 CONTINUE 

CALL MATCON 


KJ1 = 1 

00001690 

KQ1 = KJ1 + 1 

00001700 

KQ2 = KJ1 + 2 

00001710 


00001720 

DO 100 1=1, LAW 

00001730 

DO 101 J=l, LAT 

00001740 
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0095 

C 

DO 3000 IM = K J 1 , KQ2 

00001750 

00001760 

0096 


DO 3000 JM = 1, JQMAX 

00001770 

0097 


At I M , JM ) = 0. 

00001780 

0098 

3000 CONTINUE 

00001790 

0099 

C 

11=1-1 

00001800 

00001810 

0100 


12=1-2 

00001820 

0101 


Z = (FLOAT(J)— ( FLOAT ( LAT ) +1 . ) /2 • ) *H 

00001830 

0102 


NODE = LAT* I 1+ J 

00001840 

0103 


JJ1 = 3*tLAT*Il+J)-2 

00001850 

0104 


JJ2 = 3*tLAT*I2+J)-2 

00001860 

0105 


JJ3 = 3*tLAT*I2+J)-5 

00001870 

0106 


JJ4 = 3* ( L AT* I+J ) —2 

00001880 

0107 


J J5 = 3* ( LAT* 1 + J ) + 1 

00001890 

0108 


JJ6 = 3* ( LAT* I 1+ J ) + 1 

00001900 

0109 


JJ7 = 3*tLAT*I2+J)+l 

00001910 

0110 


J J 8 = 3*tLAT*Il+J)-5 

00001920 

0111 


JJ9 = 3* t L AT* I+J ) -5 

00001930 

0112 


JJ10 -3*(LAT*I1+J>— 8 

00001940 

0113 


JJ11 =3*<LAT*< I+l)+J)-2 

00001950 

0114 


JJ12 =3*(LAT*Il+J)+4 

00001960 

0115 


J J 13 =3* ( LAT* t I— 3 )+J ) —2 

00001970 

0116 

C 

JQ1 = JJ1+1 

00001980 

00001990 

0117 


JQ2 = J J 1+2 

00002000 

0118 

C 

00 102 M=l, NL AY 

00002010 

00002020 

0119 


IFtM.EQ.l. AND.J.GT . I NF ( 1 ) ) GO TO 102 

00002030 

0120 


IF(M.EQ.L) GO TO 192 

00002040 

0121 


I F C J.LE. INFtM-1) .OR.J.GT.INF(M) ) GO TO 102 

00002050 

0122 


192 IF ( J. EQ. 1 ) GO TO 200 

00002060 

0123 


IFt I.EQ.IMID.AND.J.EQ.JMID) GO TO 203 

00002070 

0124 


I F < I.EQ. IMIO+1. AND, J.EQ. JMID) GO TO 203 

00002080 

0125 


IFU.EQ.LAT> GO TO 202 

00002090 

0126 


IFt J.EQ. INF ( M) ) GO TO 201 

00002100 


C 

C 

SHOULD J EQUAL NONE OF THE ABOVE, CONTINUE ON BELOW TO STATEMENT 193 

00002110 

00002120 

0127 

c 

193 IFt I.EQ.l) GO TO 194 

00002130 

00002140 

0128 


IFU.EQ.FSW1.0R.I.EQ.FSW2) GO TO 195 

00002150 

0129 


IFt I.LT.FSW2.AND. I.GT.FSW1J GO TO 197 

00002160 

0130 


IFt I.EQ. LAW) GO TO 198 

00002170 


c 

c 

EQUILIBRIUM MATRIX TERMS FOR A SQUARE MESH, H1=H2=H3=H 

00002180 

00002190 

0131 

c 

A ( K J 1 , J J 1 > =-8.*(C66tM)+C55lM) ) 

00002200 

00002210 

0132 


A ( K J 1 , JJ2 ) = 4. *C66 ( M) 

00002220 

0133 


A t K J 1 , J J4 ) = 4 * *C66 ( M ) 

00002230 

0134 


A ( K J 1 , J J6 J = 4. *C55 ( M ) 

00002240 

0135 


A(KJ1,JJ8) = 4 ,*C55 ( M ) 

00002250 

0136 


At KJ1, JJ1+1 ) = — 8.*(C26(M) +C45 ( M ) ) 

00002260 

0137 


A(KJ1,JJ2+1) = 4 . *C26 t M ) 

00002270 

0138 


A ( K J 1 f J J4+ 1 ) = 4. *C 26 ( M > 

00002280 

0139 


AtKJlfJJ 6+ 1 ) = 4 . *C45 ( M ) 

00002290 

0140 


AtKJl, JJ8+1) = 4. *C45 ( M ) 

00002300 

0141 

c 

C = C36 1 M ) +C45 l M ) 

00002310 

00002320 
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0142 

C 

A ( K J 1 * J J3+2 ) = C 


00002330 

00002340 

0143 

A ( K J 1 1 J J5+2 ) = C 


00002350 

0144 

A ( K J 1 * J J7+2 ) = -C 


00002360 

0145 

A(KJl,JJ9+2) = -C 


00002370 

0146 

C 

X(JJ1) = 0. 


00002380 

00002390 

0147 

c 

A ( KQ1 v J J 1 ) = -8. * ( C26 ( M) +C45 { M ) ) 


00002400 

00002410 

0148 

A ( KQ1 f J J2 ) = 4 . *C26 ( M ) 


00002420 

0149 

A ( KQ1 * J J4 ) = 4.*C26tM) 


00002430 

0150 

A ( KQ1 t J J6 ) = 4.*C45(M) 


00002440 

0151 

A ( KOI » J J8 ) = 4. *C45 ( M ) 


00002450 

0152 

A ( KQ1 * J Jl+1 ) = -8 .*( C22 ( M ) +C44( M ) ) 


00002460 

0153 

A ( KQ1 y JJ2+1 ) = 4.*C22(M) 


00002470 

0154 

A (KQ1 t J J4+1 ) = 4.*C22(M) 


00002480 

0155 

A ( KQ 1 T J J 6+ 1 ) = 4.*C44(M) 


00002490 

0156 

A { KQ1 t J J8+1 ) = 4 .*C44 ( M ) 


00002500 

0157 

C 

D = C23 ( M ) +C44 ( M ) 


00002510 

00002520 

0158 

C 

A ( KQ1 » J J3+2 ) = D 


00002530 

00002540 

0159 

A ( KQ 1 f J J5+2 ) * D 


00002550 

0160 

A ( KQ1 * J J7+2 ) = -D 


00002560 

0161 

A { KQ1 1 J J9+2 ) = -D 


00002570 

0162 

C 

X ( JQ1 ) = 0. 


00002580 

00002590 

0163 

C 

A ( KQ2 » J J 3 ) = C 


00002600 

00002610 

0164 

A ( KQ2 » J J 5 ) = C 


00002620 

0165 

A(KQ2»JJ7) = -C 


00002630 

0166 

A ( KQ2 t J J9 ) = -C 


00002640 

0167 

A ( KQ2* JJ3+1 ) = D 


00002650 

0168 

A ( KQ2 t J J5+1 ) = D 


00002660 

0169 

A ( KQ2 » J J7+ 1 ) = -D 


00002670 

0170 

A ( KQ2 » J J9+ 1 ) = -D 


00002680 

0171 

A ( KQ2 1 J Jl + 2 ) = ~8.*(C44(M)+C33(M) ) 


00002690 

0172 

A ( KQ2 » J J2 + 2 ) = 4 • *C44 ( M ) 


00002700 

0173 

A ( KQ2 » J J4+2 ) = 4.*C44(M) 


00002710 

0174 

A(KQ2,JJ6+2) = 4.*C33(M) 


00002720 

0175 

A(KQ2,JJ8+2) = 4.*C33(M) 


0000273 0 

0176 

C 

XUQ2J = -4.*(C13(M)*C2 + C23<M)*DV 

+ 2.*C36(M)*C4 J*HSQ 

00002740 
0000275 0 

0177 

GO TO 102 


00002760 


C 

C FREE SURFACE MATRIX TERMS FOR 1=1 AND J 

NOT EQUAL TO 1, INF OR LAT 

00002770 
0000278 0 

0178 

C 

194 A(KJltJJl) » -3.*C66(M) 


00002790 

00002800 

0179 

A ( K J 1 f J J4 ) = 4.*C66(M) 


00002810 

0180 

A(KJ1 t JJ11) = -C661M) 


00002820 

0181 

A ( K J 1 » J J 1 + 1 ) = -3.*C26(M) 


00002830 

0182 

A ( K J 1 » J J4+1 ) = 4.4C26CM) 


00002840 

0183 

A(KJ1, JJ11+1) = -C26 ( M ) 


00002850 

0184 

A ( K J 1 » J J6+2 ) = C36 ( M ) 


00002860 

0185 

A(KJl,JJ8 + 2) = — C 36 ( M ) 


00002870 

0186 

C 

A ( KQ 1 1 J J 1 ) = — 3. J f'C26{ M) 


00002880 

00002890 

0187 

A ( KQ1 » J J4 ) =4.*C26(M) 


00002900 
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0188 



A(KQltJJll) » — C 2 6 1 M ) 

00002910 

0189 



A ( KQlt J Jl+1 ) = — 3.*C22 ( M ) 

00002920 

0190 



A(KQ1*JJ4+L) = 4.*C22(M) 

00002930 

01 * 



A(KQly JJ11+1 ) = — C22 ( M ) 

00002940 

01 



A ( KQ1 » J J6+2 ) « C231M) 

00002950 

0193 



A ( KQ1 » J J8+2 ) = — C23 ( M ) 

00002960 


C 



00002970 

0194 



A I KQ2 « J J6 ) = C45IM) 

00002980 

0? '=>5 



A ( KQ2 * J J 8 ) a — C45 ( M ) 

00002990 

0196 



A ( KQ2 * J J6+ 1 ) = C44(M) 

00003000 

0197 



A { KQ2 * J J 8+ 1 ) * — C44 ( M ) 

00003010 

0198 



A ( KQ2 * J J 1 + 2 ) = -3.*C44(M> 

00003020 

0199 



A ( KQ2 * J J4+2 ) = 4.*C44(M) 

00003030 

0200 



A ( KQ2 * J J 1 1 + 2 ) = -C44( M ) 

00003040 


C 



00003050 

0201 



CY1 = C 12 ( M ) *C3 + C22(M)*BV + C26(M)*BU 

00003060 

0202 



CY2 = C12(M)*C2 + C22 ( M ) *DV + 2.*C26(M)*C4 

00003070 

0203 



CXY 1 = C 16 ( M ) *C3 + C26(M)*BV + C66<M)*BU 

00003080 

0204 



CXY2 = C16(M)*C2 + C26(M)*0V + 2 .*C66 ( M ) *C4 

00003090 


C 



00003100 

0205 



X(JJ1) = —2 • *H* ( CXY1 + CXY2*Z) 

00003110 

0206 



X(JQ1) = -2.*H*(CY1 + CY2*Z) 

00003120 

0207 



X ( JQ2 ) = 0. 

00003130 

0208 



GO TO 102 

00003140 


C 



00003150 

0209 


195 

HI = H 

00003160 

0210 



H2 = FLOAT t K ) *H 

00003170 

0211 



H3 = H 

00003180 


c 



00003190 

0212 



IFd.NE.FSW2) GO TO 196 

00003200 

0213 



HI = FLOAT(K)*H 

00003210 

0214 



H2 = H 

00003220 


c 



00003230 

0215 


196 

CONTINUE 

00003240 

0216 



HH = H2/H1 

00003250 

0217 



HR = HH/d.+HH) 

00003260 

0218 



HH1 = H1/H3 

00003270 

0219 



HH2 = H2/H3 

00003280 

0220 



HH3 = H1*H2 

00003290 

0221 



HMU = HH1*HH2 

00003300 

0222 



GO TO 199 

00003310 


c 



00003320 

0223 


197 

HI = FLOAT ( K > *H 

00003330 

0224 



H2 = HI 

00003340 

0225 



H3 = H 

00003350 

0226 



GO TO 196 

00003360 


c 



00003370 


c 

FREE SURFACE MATRIX TERMS FOR I=LAW AND J NOT EQUAL TO 1, 

INF OR LAT 00003380 


c 



00003390 

0227 


198 

A(KJltJJl) = 3.*C66(M) 

00003400 

0228 



A ( K J 1 » J J2 ) = -4.*C66{M) 

00003410 

0229 



A ( K J 1 1 J J 1 3 ) a C66 ( M ) 

00003420 

0230 



A ( K J lyJJl+1) = 3.*C26(M) 

00003430 

0231 



A ( K J 1 » J J2+ 1 ) = -4.*C26(M) 

00003440 

0232 



A ( K J 1 » J J 1 3+ 1 ) * C26 ( M ) 

00003450 

0233 



A ( K J 1 * J J6+2 ) = C36 ( M ) 

00003460 

0234 



A { K J 1 * J J8+2 ) = -C36(M) 

00003470 


C 00003480 
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0235 


A ( KQ1 , JJ 1 ) = 

3. *C26 ( M ) 

00003490 

0236 


A ( KQ1 y J J2 ) = 

— 4.*C26 C M ) 

00003500 

0237 


A { KQ1 * JJ 13 ) 

= C26 { M ) 

00003510 

0238 


A ( KQ1 , J J 1+1 ) 

= 3.*C22(M) 

00003520 

0239 


A(KQ1, JJ2+1) 

= — 4. *C22 ( M) 

00003530 

0240 


A ( KQ1 , J J13+1 

) = C221M) 

00003540 

0241 


A ( KQ1 » JJ6+2 ) 

= C23 ( M ) 

00003550 

0242 

C 

AiKQl, JJ8+2) 

= -C23(M) 

00003560 

00003570 

0243 


A ( KQ2 » J J6 ) = 

C45 ( M ) 

00003580 

0244 


A ( KQ2 » J J8 ) = 

— C45 ( M ) 

00003590 

0245 


A ( KQ2 * J J6+1 ) 

= C44 ( M ) 

00003600 

0246 


A(KQ2, JJ8+1) 

= -C44 ( M ) 

00003610 

0247 


A(KQ2,JJl+2) 

= 3 »*C44 { M ) 

00003620 

0248 


A(KQ2, JJ2+2) 

= — 4.*C44 ( M ) 

00003630 

0249 


A ( KQ2 « J J 1 3+2 ) = C44 ( M ) 

00003640 


C 



00003650 

0250 


CY1 = C 1 2 ( M ) *C3 + C22 ( M ) * BV + C26(M)*BU 

00003660 

0251 


CY2 = C 12 ( M ) *C2 + C22 ( M ) *DV + 2.*C26<M)*C4 

00003670 

0252 


CXY1 = C16<M)*C3 + C26(M)*BV + C66(M)*BU 

00003680 

0253 


CXY2 = C16(M)*C2 + C26(M)*DV + 2.*C66{ M) *C4 

00003690 


c 



00003700 

0254 


X(JJ1) = -2. 

*H*<CXY1 + CXY2*2 ) 

00003710 

0255 


X(J01) = -2. 

{ C Y 1 + CY2*Z) 

00003720 

0256 


X { JQ2 ) = 0. 


00003730 

0257 

c 

GO TO 102 


00003740 

00003750 


c 

EQUILIBRIUM MATR 

IX TERMS FOR A VARIABLE MESH, HI, H2 

, H3 INDEPENDENT 00003760 


c 



00003770 

0258 


199 A ( KJ 1 , J J 1 ) = 

-2.*(C66(M)+HMU*C55(M) ) 

00003780 

0259 


A ( K J 1 » J J2 ) = 

2.*HR*C66( M) 

00003790 

0260 


A (K J1 , JJ4) = 

2.*C66(M)/( l.+HH) 

00003800 

0261 


A ( K J 1 , J J6 ) = 

HMU*C55 ( M ) 

00003810 

0262 


A{KJ1,JJ8) = 

HMU*C 55 { M ) 

00003820 

0263 


A { K Jl, JJ1+1 ) 

= -2.*(C26(M) + HMU*C45< M) ) 

00003830 

0264 


A ( KJ 1 » J J2+ 1 ) 

- 2.*HR*C26(M) 

00003840 

0265 


A(KJ1, JJ4+1) 

= 2.*C26<M)/( l.+HH) 

00003850 

0266 


A (KJ1 , JJ6+1 ) 

= HMU*C45(M) 

00003860 

0267 

c 

A(KJ1, JJ8+1) 

= HMU*C45{M) 

00003870 

00003880 

0268 


C - HH 1*HR* ( C36 ( M ) +C45 ( M) ) /2. 

00003890 


c 



00003900 

0269 


A ( K J 1 t J J3+2 ) 

= C 

00003910 

0270 


A { KJ 1 , J J5+2 ) 

= c 

00003920 

0271 


AtK'Jl, JJ7+2) 

= -c 

00003930 

0272 

c 

A ( K J 1 » J J9+2 ) 

= -c 

00003940 

00003950 

0273 


A(KQ1 t JJ1) - 

— 2 ( C26 ( M ) +HMU*C45 ( M ) ) 

00003960 

0274 


A ( KQ1 , J J2 ) = 

2.*HR*C26( M) 

00003970 

0275 


A(KQ1 , JJ4) = 

2.*C26(M}/{ l.+HH) 

00003980 

0276 


A t KQ1 , J J6 ) = 

HMU*C45< M) 

00003990 

0277 


A ( KQ1 f J J 8 ) = 

HMU*C45 ( M ) 

00004000 

0278 


A ( KQ1, JJ1+1 ) 

= -2.*(C22(M)+HMU*C44<M) ) 

00004010 

0279 


A ( KQ1 * J J2+ 1 ) 

= 2.*HR*C22(M) 

00004020 

0280 


A ( KQlt J J4+1 ) 

= 2.*C22(M)/( l.+HH) 

00004030 

0281 


A { KQ1 ,JJ6+1) 

= HMU*C44(M) 

00004040 

0282 


A(KQ1, JJ8+1) 

= HMU*C44(M) 

00004050 


C 00004060 
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0283 


D = HH1*HR*(C23(M)+C44(M) ) /2. 

OC 004070 


C 


00004080 

0284 


A ( KQ1 * J J3+2 ) = D 

00004090 

0285 


A ( KQ1 * J J5+2 ) = D 

00004100 

0286 


A ( KQ1 » J J7+2 ) * -D 

00004110 

0287 


A ( KQ1 * J J9+2 ) = -D 

00004120 


C 


00004130 

0288 


A ( KQ2 » J J3 ) = C 

00004140 

0289 


A ( KQ2 ♦ J J 5 ) = C 

00004150 

0290 


A ( KQ2* JJ7 ) = -C 

00004160 

0291 


A ( KQ2 * J J9 ) = -C 

00004170 

0292 


A ( KQ2 1 J J3+1 1 = D 

00004180 

0293 


A ( KQ2 » J J5+ 1 ) = D 

00004190 

0294 


A ( KQ2 » J J7+ 1 ) * -D 

00004200 

0295 


A ( KQ2 * J J9+1 ) = -D 

00004210 

0296 


A ( KQ2» J J 1+2 ) * -2.*(C44(M)+HMU*C33(M) ) 

00004220 

0297 


A ( KQ2 t J J2+2 ) * 2.*HR*C44(M) 

00004230 

0298 


A ( KQ2* J J4+2 ) * 2.*C44(M)/( l.+HH) 

00004240 

0299 


A ( KQ2 r J J 6+2 ) = HMU*C33(M) 

00004250 

0300 


A ( KQ2* J J8+2 ) = HMU*C33(M) 

00004260 


C 


00004270 

0301 


XIJJ1) * 0. 

00004280 

0302 


X { JQ1 ) = 0. 

00004290 

0303 


X ( JQ2 ) = — HH3* (C13(M)*C2 + C23(M)*DV + 2 . *C36 ( M ) *C4 ) 

00004300 

0304 


GO TO 102 

00004310 


C 


00004320 

0305 


200 I F ( I . EO. 1 ) GO TO 210 

00004330 

0306 


I F ( I.EQ.LAW) GO TO 211 

00004340 


C 


00004350 


c 

FREE SURFACE MATRIX TERMS FOR I BETWEEN 1 AND LAW AND J=l. 

00004360 


c 


00004370 

0307 


A { K J 1 » J J 1 ) = -3.*C55(M) 

00004380 

0308 


A ( K J 1 r J J6 > = 4.*C55(M) 

00004390 

0309 


A ( K J 1 » J J12 ) = -C55 ( M ) 

00004400 


c 


00004410 

0310 


A ( K J 1 « JJ 1 + 1 ) = -3.*C45(M) 

00004420 

0311 


A ( K J 1 1 J J 6+ 1 ) = 4.*C45(M) 

00004430 

0312 


A(KJ1, JJ12+1) = -C45(M) 

00004440 


c 


00004450 

0313 


A(KQ1,JJ1) = -3„*C45(M> 

00004460 

0314 


A ( KQ1 » J J6 ) = 4 .*C45 ( M ) 

00004470 

0315 


A(K01»JJ12) = -C45 ( M ) 

00004480 


c 


00004490 

0316 


A(KQlyJJl+l> = — 3. *C44 ( M ) 

00004500 

0317 


A 1 KQ1 y J J6+ 1 ) = 4.*C44(M) 

00004510 

0318 


A ( KQly JJ12+1 ) = — C44 ( M ) 

00004520 


c 


00004530 

0319 


A ( KQ2 * J J 1 + 2 ) = -3.*C33(M) 

00004540 

0320 


A ( KQ2 * J J6+2 ) = 4.*C33(M) 

00004550 

0321 


AIKQ2, JJ12+2) = -C33 ( M ) 

00004560 


c 


00004570 

0322 


CZ1 = C13(M)*C3 + C23{ M) *BV + C36(M)*BU 

00004580 

0323 


CZ2 = C13(M)*C2 + C23(M)*DV + 2.*C36(M)*C4 

00004590 


c 


00004600 

0324 


X( JJ1 ) = 0. 

00004610 

0325 


X ( J01 ) = 0. 

00004620 

0326 


X ( JQ2 ) = -2.*H*(CZ1 + CZ2*Z ) 

00004630 


c 


00004640 
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0327 



I F C I.EQ.FSW1 ) GO TO 206 


00004650 

0328 



I F C I.EQ.FSW2) GO TO 206 


00004660 

0329 



IF( I.GT.FSW1. AND. I.LT.FSW2) GO TO 

209 

00004670 


C 




00004680 


c 

IF 

I IS BETWEEN 1 AND FSWi OR BETWEEN 

FSW2 AND LAW, CONTINUE BELOW 

00004690 


c 




00004700 

0330 



A( KJ1, JJ2+2) = -C45IMJ 


00004710 

0331 



A ( K J 1 * J J4+2 ) = C45 ( M ) 


00004720 


c 




00004730 

0332 



A ( KQ1 » J J2 + 2 ) = — C44( M.) 


00004740 

0333 



A<KQl,JJ4+2) = C44 ( M ) 


00004750 


c 




00004760 

0334 



A ( KQ2 , JJ2 ) = — C36 { M ) 


00004770 

0335 



A ( KQ2 * J J4 ) = C36 ( M ) 


00004780 

0336 



A ( KQ2, J J2+ 1 ) = — C23 ( M ) 


00004790 

0337 



A ( KQ2 * J J4+1 ) = C23 ( M ) 


00004800 

0338 



GO TO 102 


00004810 


c 




00004820 


c 

CASE WHERE I=FSW1 OR FSW2 AND J=1 


00004830 


c 




00004840 

0339 


206 

XK a FLOAT ( K ) 


00004850 

0340 



D1 = 2 • * ( XK— 1 . ) /XK 


00004860 

0341 



D2 = 2.*XK/( XK+1. ) 


00004870 

0342 



D3 = 2 • / ( {XK+1. ) *XK) 


00004880 


c 




00004890 

0343 



IF ( I .EQ.FSW2 ) GO TO 207 


00004900 


c 




00004910 

0344 



A(KJltJJl+2) = D1*C45(M) 


00004920 

0345 



A ( K J 1 , J J2 + 2 ) a -D2*C45(M) 


00004930 

0346 



A { K J 1 » J J4+2 ) = D3*C45<M) 


00004940 


c 




00004950 

0347 



A { KOI , J J 1 + 2 ) = D 1*C44 ( M ) 


00004960 

0348 



A { KQ1 * JJ2 + 2 ) = -D2*C44(M) 


00004970 

0349 



A ( KQ1 » J J4+2 ) = D3*C44(M) 


00004980 


c 




00004990 

0350 



A ( KQ2 * J J 1 ) = D1*C36(M) 


00005000 

0351 



A ( K02 1 JJ2 ) = -D2*C 36 { M ) 


00005010 

0352 



A ( KQ2 t J J4 ) a D3*C36 { M ) 


00005020 


c 




00005030 

0353 



A(KQ2,JJ1+1) a D1*C23(M) 


00005040 

0354 



A(KQ2,JJ2+1) = -D2*C23 ( M ) 


00005050 

0355 



AIKQ2,JJ4+1) = D3*C23(M) 


00005060 

0356 



GO TO 102 


00005070 


c 




00005080 

0357 


207 

A(KJl,JJl+2) = -D1*C45(M) 


00005090 

0358 



A ( K J 1 » J J2 + 2 ) = -D3*C45(M) 


00005100 

0359 



A ( K J 1 v J J4+2 ) = D2*C45(M) 


00005110 


c 




00005120 

0360 



A ( KOI » J J 1 + 2 ) = -D1*C44(M> 


00005130 

0361 



A(KQl t JJ2+2) = — D3*C44 ( M ) 


00005140 

0362 



A(KQl,JJ4+2) = D2*C44(M) 


00005150 


c 




00005160 

0363 



A ( K02 » J J 1 ) = -D1*C36(M) 


00005170 

0364 



A(KQ2,JJ2) = -D3^C36(M) 


00005180 

0365 



A ( KQ2 y J J4 ) = D2*C36(M) 


00005190 


c 




00005200 

0366 



A(KQ2 t JJl+l) = -D1*C23(M) 


00005210 

0367 



A(KQ2tJJ2+1) = -D3«C23(M) 


00005220 
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0368 



A ( KQ2, J J4+1 ) = D2*C23(M> 

00005230 

0369 



GO TO 102 

00005240 


C 



00005250 


c 

CASE WHERE I IS BETWEEN FSW1 AND FSW2 ANO J=1 

00005260 


c 



00005270 

0370 


209 

XK = FLOAT ( K ) 

00005280 

0371 



A ( K J 1 » J J2+2 ) - -C45(M)/XK 

00005290 

0372 



A ( K J 1 , J J4+2 ) = C45 ( M ) / XK 

00005300 


c 



00005310 

0373 



AtKQl, JJ2+2) * — C44 { M } /XK 

00005320 

0374 



AtKQl , J J4+2 ) = C441MJ/XK 

00005330 


c 



00005340 

0375 



A ( KQ2 » J J2 ) = — C36 { M ) /XK 

00005350 

0376 



A ( KQ2* J J4 ) = C 36 ( M } /XK 

00005360 


c 



00005370 

0377 



A I KQ2 , J J2+ 1 ) = — C2 3 CM) /XK 

00005380 

0378 



A ( KQ2 » J J4+ 1 ) = C23<M)/XK 

00005390 

0379 



GO TO 102 

00005400 


c 



00005410 


c 

FREE SURFACE MATRIX TERMS FOR I=J=1 

00005420 


c 



00005430 

0380 


210 

A(KJltJJl) = -3.*C66(M) 

00005440 

0381 



A ( K J 1 * J J4 ) = 4.*C66(M) 

000Q5450 

0382 



A(KJltJJll) = -C 66 (M) 

00005460 

0383 



A(KJ1 t JJ1+1) = -3.*C26(M) 

00005470 

0384 



A ( K J 1 y J J4+ l ) = 4.*C26<M> 

00005480 

0385 



A t K J l 9 J J 1 1+ 1 ) = -C26(M) 

00005490 

0386 



A ( K J 1 9 J J 1+2 ) = -3.*C36<M) 

00005500 

0387 



A (KJ1 9 JJ6+2 ) * 4.*C361M) 

00005510 

0388 



A ( K J 1 9 J J 12+2 ) = -C36IM) 

00005520 


c 



00005530 

0389 



A { KQ 1 9 J J 1 ) = -3.*C26(M) 

00005540 

0390 



A ( KOI 9 J J4 ) = 4.*C26(M) 

00005550 

0391 



A(KOlfJJll) = -C261M) 

00005560 

0392 



A ( KQ1 9 J J 1+ 1 ) = -3.*C22tM) 

00005570 

0393 



AtKQl, J J4+ 1 ) = 4.*C22(M) 

00005580 

0394 



AtKQl, JJ11+1) = — C 22 ( M ) 

00005590 

0395 



A ( KQ 1 1 J J 1 + 2 ) = -3.*C23<M) 

00005600 

0396 



AtKQl, JJ6+2) = 4.*C23(M> 

00005610 

0397 



AtKQl, JJ12 + 2) = -C 2 3 ( M ) 

00005620 


c 



00005630 

0398 



A(KQ2,JJ1) = — 3« *C45 ( M) 

00005640 

0399 



A ( KQ2 , J J 6 ) = 4 • *C45 ( M ) 

00005650 

0400 



A ( KQ2 , J J 1 2 ) = — C45 ( M ) 

00005660 

0401 



A ( KQ2 , J Jl + 1 ) = -3.*C44(M) 

00005670 

0402 



A t KQ2 , J J 6 + 1 ) = 4.*C44(M) 

00005680 

0403 



A(KQ2,JJ12+1 ) = —C44 ( M ) 

00005690 

0404 



A { KQ2 , J J 1+2 ) = -3.*C44(M) 

00005700 

0405 



A ( KQ2 , J J4+2 ) = 4»*C44 ( M ) 

00005710 

0406 



A ( KQ2 , J J 1 1+2 ) = -C44{ M ) 

00005720 


c 



00005730 

0407 



CY1 = C12(M)*C3 + C22tM)*BV + C26(M)*BU 

00005740 

0408 



CY2 = C12tM)*C2 + C22tM)*DV + 2 . *C26 t M ) *C4 

00005750 

0409 



CXY 1 = C16(M)*C3 + C26(M)*BV + C 66 (M)*BU 

00005760 

0410 



CXY2 = C 1 6 ( M ) *C2 + C26(M)*DV + 2 .*C 66 l M ) *C4 

00005770 


c 



00005780 

0411 



X(JJl) = -2.*H*tCXYl + CX Y2*Z ) 

00005790 

0412 



X(JQl) = -2.*H*(CY1 + CY2*Z) 

00005800 
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0413 



X ( JQ2 ) * 0. 

00005810 

0414 



GO TO 102 

00005820 


C 



00005830 


c 

FREE SURFACE MATRIX TERMS FDR J=1 AND I=LAW 

00005-840 


c 



00005850 

0415 


211 

A(KJltJJl) = 3.*C66(M) 

00005860 

0416 



A ( K J 1 f J J2 ) = -4.*C66(M) 

00005870 

0417 



A ( K J 1 t J J 13 ) = C66 ( M ) 

00005880 

0418 



A ( KJ 1 1 J Jl+1 ) = 3.*C26(M) 

00005890 

0419 



A(KJ1 » J J2+1 ) = -4.*C26(M) 

00005900 

0420 



A ( K J 1 1 J J 1 3+ 1 ) = C26 ( M ) 

00005910 

0421 



A ( KJ 1 f J J 1+2 ) = -3.*C36 { M ) 

00005920 

0422 



A(KJl,JJ6+2) = 4.*C36<M) 

00005930 

0423 



A < K J 1 * J J 12+2 ) * — C36 ( M ) 

00005940 


c 



00005950 

0424 



A ( KQ1 v J J 1 ) = 3.*C26<M> 

00005960 

0425 



A { KQ 1 1 J J2 ) = -4.*C26{ M) 

00005970 

0426 



A ( KQ1 t J J 1 3 ) = C26 ( M ) 

00005980 

0427 



A ( KQ1 1 J J 1+ 1 ) = 3. *C22 ( M ) 

00005990 

0428 



A { KQ1 r J J2+ 1 ) = — 4.*C22 ( M ) 

00006000 

0429 



A { KQ1 , J J13+1 ) = C22 ( M) 

00006010 

0430 



A { KQ1 1 J J 1 + 2 ) = -3.*C23(M) 

00006020 

0431 



A ( KQ1 f J J6+2 ) = 4.*C23(M) 

00006030 

0432 



A(KQl t JJ12+2) = -C23 ( M ) 

00006040 


c 



00006050 

0433 



A ( KQ2 t J J 1 ) = —3 • *C45 ( M ) 

00006060 

0434 



A ( KQ2 » J J6 ) = 4.*C45(M) 

00006070 

0435 



A ( KQ2 f J J 12 ) = — C45 ( M ) 

00006080 

0436 



A ( KQ2 * J J 1+ 1 ) * -3. *C44 ( M ) 

00006090 

0437 



A ( KQ2 * J J6+ 1 ) = 4.*C44(M) 

00006100 

0438 



A ( KQ2» J J12+ 1 ) = -C44 ( M ) 

00006110 

0439 



A { KQ2 t J J 1 + 2 ) = 3.*C44(M) 

00006120 

0440 



A ( KQ2 * J J2+2 ) » -4.*C44(M) 

00006130 

0441 



A ( KQ2 » J J 13 + 2 ) - C44 ( M ) 

00006140 


c 



00006150 

0442 



CY1 = C 12 ( M ) *C3 + C22(M)*BV + C26(M)*BU 

00006160 

0443 



CY2 a C12(M)*C2 + C22 ( M ) *DV + 2.*C26(M)*C4 

00006170 

0444 



CXY1 = C 1 6 ( M ) *C3 + C26(M)*BV + C66(M)*BU 

00006180 

0445 



CXY2 = C16<M)*C2 + C26(M)*DV + 2.*C66(M)*C4 

00006190 


c 



00006200 

0446 



X(JJ1) = — 2. #H* ( CX Y 1 + CXY2*Z) 

00006210 

0447 



X { JQ1 ) = -2.*H*(CY1 + CY2*Z) 

00006220 

0448 



X ( JQ2 ) = 0. 

00006230 

0449 



GO TO 102 

00006240 


c 



00006250 

0450 


201 

P = M+l 

00006260 

0451 



I F { I.EQ.l) GO TO 220 

00006270 

0452 



I F { I.EQ.FSW1) GO TO 221 

00006280 

0453 



IF( I.LT.FSW2. AND. I.GT.FSW1 ) GO TO 222 

00006290 

0454 



IFU.EQ.FSW2) GO TO 221 

00006300 

0455 



I F ( I.EQ.LAW) GO TO 223 

00006310 


c 



00006320 


c 

MATRIX TERMS AT INTERFACE FOR I BETWEEN 1 AND FSW1 

. OR FSW2 AND LAW 00006330 


c 



00006340 

0456 



A(KJlfJJl) * 3.*(C55(M)+C55(P) ) 

00006350 

0457 



A(KJl r JJ6) = -4.*C55(P) 

00006360 

0458 



A { K J 1 1 J. J 8 ) = -4.*C55(M) 

00006370 

0459 



A(KJl f JJ10) = C55 ( M) 

00006380 
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0460 



A(KJ1,JJ12) = C55 ( P ) 

00006390 


C 



00006400 

0461 



A ( K J 1 , J J 1+ 1 ) = 3.*<C45(M)+C45(P> > 

00006410 

0462 



A ( K J 1 • J J6+1 ) = -4.*C45(P> 

00006420 

0463 



A ( K J 1 , J J8+1 ) = -4.*C45(M> 

00006430 

0464 



AtKJl, JJ10+1) = C45 < M ) 

00006440 

0465 



A ( KJ 1 » J J 12+1 ) = C45 ( P ) 

00006450 


C 



00006460 

0466 



A ( KJ 1 , J J2+2 ) * C45 ( P ) — C45 ( M ) 

00006470 

0467 



A ( KJ 1 , J J4+2 ) = C45tM)-C45< P) 

00006480 


c 



00006490 

0468 



A(KQ1»JJ1) * 3.*(C45<M)+C45(P) ) 

00006500 

0469 



A ( KQ1 , J J6 ) = -4.*C45(P) 

00006510 

0470 



AIKQ1.JJ8) » -4.*C45(M) 

00006520 

0471 



A ( KQ1 « J J 10 ) * C45 ( M ) 

00006530 

0472 



A ( KOI t J J12 ) = C45 ( P ) 

00006540 


c 



00006550 

0473 



A ( KOI , J J 1+1 ) = 3, * { C44 ( M ) +C44( P ) ) 

00006560 

0474 



A ( KQ1 ♦ J J6+1 ) = -4.*C44(P) 

00006570 

0475 



A { KQ1 , J J8+ 1 ) = -4* #C44( M ) 

00006580 

0476 



A ( KQ1 , J J 10+1 ) = C44( M ) 

00006590 

0477 



A(KQ1, JJ12+1) = C44{P) 

00006600 


c 



00006610 

0478 



A ( KQ1 1 J J2+2 ) = C44{ P J -C44 ( M ) 

00006620 

0479 



A { KQ1 , J J4+2 ) = C44 { M ) — C44 ( P ) 

00006630 


c 



00006640 

0480 



A { KQ2 * J J2 ) = C36 I P ) — C36 ( M ) 

00006650 

0481 



A ( KQ2 » J J4 ) = C36 ( M ) -C36 ( P ) 

00006660 


0 



00006670 

0482 



A ( KQ2 * J J2+ 1 ) = C23 ( P ) — C23 ( M ) 

00006680 

0483 



A(KQ2*JJ4+1) = C23(M)-C23(P) 

00006690 


c 



00006700 

0484 



A ( KQ2 * J J 1 + 2 > = 3.*tC33(M)+C33(P) ) 

00006710 

0485 



A ( KQ2* JJ6+2 ) = ~4.*C33(P) 

00006720 

0486 



A(KQ2,JJ8+2) = ~4 ■ *C3 3 ( M ) 

00006730 

0487 



A(KQ2, JJ10+2) = C33 ( M) 

00006740 

0488 



A ( KQ2 f J J 1 2 + 2 ) = C33 ( P ) 

00006750 


c 



00006760 

0489 



CZ1 = < C 13 < P ) — C 13 < M ) )*C3 + (C23(P)-C23(M) )*BV + ( C36 ( P ) -C36 ( M )) *BU00006770 

0490 



CZ2 = ( C 1 3 { P ) — C 1 3 ( M ) ) *C2+ ( C23 ( P ) -C 23 ( M ) ) *0V+2 . * { C36 ( P ) -C36 ( M ) ) *C4 

00006780 


c 



00006790 

0491 



X{ JJ1) a o. 

00006800 

0492 



XIJQ1) = 0. 

00006810 

0493 



X ( J02 ) a 2.*H*<CZ1 + CZ2*Z ) 

00006820 

0494 



GO TO 102 

00006830 


c 



00006840 


c 

FREE 

1 SURFACE MATRIX TERMS AT ANY INTERFACE WHERE 1=1 AND J=INF OR AT 

00006850 


c 

THE 

FREE SURFACE POINT 1=1, J= LAT 

00006860 


c 



00006870 

0495 


220 

A(KJ1,JJ1) = -3,*C66(M) 

00006880 

0496 



A(KJ1,JJ4) = 4.*C66(M) 

00006890 

0497 



A ( K J 1 , J J 1 1 ) = -C66 ( M) 

00006900 


c 



00006910 

0498 



A ( K J 1 , J J 1+1 ) = -3.*C26(M) 

00006920 

0499 



A ( K J 1 , J J4+ 1 ) = 4.^C26(M) 

00006930 

0500 



A(KJ1, JJ11+1) = -C26(M) 

00006940 


c 



00006950 

0501 



A ( K J 1 , J J 1+2 ) = 3.*C36(M) 

00006960 
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0502 


A < K J 1 * J J8+2 ) = -4.*C36(M) 

00006970 

0503 


A ( K J 1 1 J J 10+2 ) = C361M) 

00006980 


C 


00006990 

0504 


A ( KQ1 » J J 1 ) = -3.*C26(M) 

00007000 

0505 


A ( KQ1 y JJ4) =4.*C26(M) 

00007010 

0506 


A (KQ1 v J J1 1 ) = — C26 ( M ) 

00007020 


C 


00007030 

0507 


A ( KOI » J J 1+1 ) = -3.*C22(M) 

00007040 

0508 


A { KQlt JJ4+1 ) = 4.*C22(M> 

00007050 

0509 


A ( KQ1 * J J 1 1+1 ) = — C22 ( M ) 

00007060 


c 


00007070 

0510 


A ( KQ1 » J J 1+2 ) = 3.*C23(M) 

00007080 

0511 


A ( KQl * J J8+2 ) = -4.*C23(M) 

00007090 

0512 


A { KQ1 , J J 10+2 ) = C23 ( M ) 

00007100 


c 


00007110 

0513 


A ( KQ2 « J J 1 ) = 3.*C45(M) 

00007120 

0514 


A ( KQ2 » J J 8 ) = —4* *C45 ( M ) 

00007130 

0515 


A ( KQ2 f J J10 ) = C45 ( M ) 

00007140 


c 


00007150 

0516 


A(KQ2,JJ1+1> * 3.*C44(M> 

00007160 

0517 


A ( K02 1 J J8+ 1 ) = -4.*C44(M) 

00007170 

0518 


A (K02 * J J10+1 ) = C44 ( M ) 

00007180 


c 


00007190 

0519 


A ( KQ2 « J J 1 + 2 ) = —3 .*C44 ( M ) 

00007200 

0520 


A { KQ2 ♦ J J4+2 ) = 4.*C44(M) 

00007210 

0521 


A ( KQ2 t J J 11 + 2) = —C44 ( M ) 

00007220 


c 


00007230 

0522 


CY1 = C 1 2 ( M ) *C3 + C22 ( M ) =*BV -*• C26<M)*BU 

00007240 

0523 


CY2 = C 12 ( M ) *C2 + C 22 ( M ) *DV + 2.*C26< M ) *C4 

00007250 

0524 


CXY1 = C 1 6 { M ) *C3 + C26 { M ) *BV + C66(M)*BU 

00007260 

0525 


CXY2 = C16(M)*C2 + C26{M)*DV + 2.*C66(M)*C4 

00007270 


c 


00007280 

0526 


X(JJ1) = -2 • ( C X Y 1 + CXY2*Z) 

00007290 

0527 


X(JQ1) = -2.*H*(CY1 + CY2*Z) 

00007300 

0528 


X ( JQ2 ) = 0, 

00007310 

0529 


GO TO 102 

00007320 


c 


00007330 


c 

MATRIX TERMS AT THE INTERFACE FOR J=INF AND I=FSW1 OR 

I = F SW 2 00007340 


c 


00007350 

0530 


221 XK = FLOAT ( K ) 

00007360 

0531 


D1 = { XK— 1 • ) /XK 

00007370 

0532 


D2 = XK / ( XK+ 1 • ) 

00007380 

0533 


D3 = 1./ ( (XK+1. )*XK) 

00007390 


c 


00007400 

0534 


A(KJlyJJl) = 3.*<C55(M)+C55( P) ) 

00007410 

0535 


A(KJ1,JJ6> * —4 . *C55 ( P ) 

00007420 

0536 


A(KJlfJJ8) = -4.*C55(M) 

00007430 

0537 


A ( K J 1 f J J 10 ) = C55 ( M ) 

00007440 

0538 


A ( K J 1 f J J 12 ) = C55 ( P ) 

00007450 


c 


00007460 

0539 


A ( K J 1 » J J 1+ 1 ) = 3.*< C45< M)+C45( P) > 

00007470 

0540 


A ( K J 1 1 J J 6+ 1 ) = -4.*C45(P> 

00007480 

0541 


A(KJ1»JJ8+1) = -4.*C45(M) 

00007490 

0542 


A ( K J 1 * J J 10+ 1 ) = C45 ( M ) 

00007500 

0543 


A(KJ1, JJ12+1) = C45 ( P ) 

00007510 


c 


00007520 

0544 


A(KQltJJl) = 3. * { C45 ( M ) + C45 ( P ) ) 

00007530 

0545 


A ( KQ 1 * J J 6 ) = -4„*C45(P) 

00007540 
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0546 


A { KQ1 * J J8 ) = -4.*C45(M> 

00007550 

0547 


A t K.Q1 * J J 10 ) = C45 { M ) 

00007560 

0548 


A ( KQ1 v J J 12 ) = C45 ( P ) 

00007570 


C 


00007580 

0549 


A ( KQ 1 f J J 1+ 1 ) = 3 •* ( C44 ( M ) +C44 ( P ) 1 

00007590 

0550 


A ( KQ 1 t J J6+1 ) = -4.*C44(P) 

00007600 

0551 


A (KQ1 » JJ8+1 ) = — 4 . *C44 ( M ) 

00007610 

0552 


A ( KQ1 * JJ10+1 ) = C44 ( M ) 

00007620 

0553 


A(KQ1, JJ12+1) = C44IP) 

00007630 


C 


00007640 

0554 


A (KQ2 f J Jl + 2 ) = 3.*(C33(M)+C33(P) ) 

00007650 

0555 


A { KQ2 » J J6+2 ) - -4.*C33(P> 

00007660 

0556 


A { KQ2 » J J8+2 ) = -4.*C33(M) 

00007670 

0557 


A ( KQ2» JJ10+2 ) = C 33 ( M ) 

00007680 

0558 


MKQ2,JJ12+21 = C33IPI 

00007690 


c 


00007700 

0559 


C21 = (C13(P)-C13(M) ) *C 3 + ( C23 < P ) -C23 I M ) ) *BV 

+ ( C36< P ) -C36( M ) ) *BU00007710 

0560 


CZ2 = (C13I PJ-C13CM) )*C2+(C23I P)-C23(M) )*DV+2. 

*(C36(P)-C36<M) )*C4 00007720 


c 


00007730 

0561 


X(JJ1) = 0, 

00007740 

0562 


X { JQ1 ) = 0. 

00007750 

0563 


X ( JQ2) = 2.*H*(CZ1 + CZ2*Z ) 

00007760 


c 


00007770 

0564 


C=C45(M)-C45(P) 

00007780 

0565 


D=C44(M)-C44( P) 

00007790 

0566 


E=C23(M)-C23(P) 

00007800 

0567 


CC-C36(M)-C36(P) 

00007810 


c 


00007820 

0568 


I F ( I.E0.FSW2) GO TO 227 

00007830 


c 


00007840 

0569 


A ( K J 1 f J J 1 + 2 ) = 2.*Dl*C 

00007850 

0570 


A(KJl,JJ2+2) * -2.*D2*C 

00007860 

0571 


A ( K J 1 1 J J4+2 ) - 2.*D3*C 

00007870 


c 


00007880 

0572 


AtKQl, JJ1+2) = 2.*D1*D 

00007890 

0573 


A { KOI t J J2 + 2 ) = -2.*D2*D 

00007900 

0574 


A ( KQ1 ♦ J J4+2 ) = 2.*D3*D 

00007910 


c 


00007920 

0575 


A { K02 y J J 1 ) = 2 .*D1*CC 

00007930 

0576 


A ( K02 » J J2 ) = -2.*D2*CC 

00007940 

0577 


A ( KQ2 * J J4 ) = 2.*D3*CC 

00007950 


c 


00007960 

0578 


A(KQ2»JJ1+1) * 2.*Dl*E 

00007970 

0579 


A 1 KQ2 t J J2+ 1 ) = -2.*D2*E 

00007980 

0580 


A 1 KQ2 » J J4+ 1 ) = 2.*D3*E 

00007990 

0581 


GO TO 102 

00008000 


c 


00008010 

0582 


227 A(KJl,JJl+2) - -2 « *Dl*C 

00008020 

0583 


A ( K J 1 » J J 2 + 2 ) = -2.*D3*C 

00008030 

0584 


A ( K J 1 » J J4+2 ) = 2.*D2*C 

00008040 


c 


00008050 

0585 


A(KQ1,JJ1 + 2J = -2.*D1*D 

00008060 

0586 


A ( KQ 1 » J J2 + 2 ) = -2.*D3*D 

00008070 

0587 


A(KQl»JJ4+2) = 2.*D2*D 

00008080 


c 


00008090 

0588 


A ( KQ2 » J J 1 ) = ^2.*D1*CC 

00008100 

0589 


A ( KQ2 y J J2 ) = -2.*D3*CC 

00008110 

0590 


A ( KQ2 t J J4 ) = 2.*D2*CC 

00008120 
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C 



00008130 

0591 


A ( KQ2 * J J 1+1 ) = -2.*D1*E 


00008140 

0592 


A ( KQ2 » J J2+ 1 ) = -2.*D3*E 


00008150 

0593 


A ( KQ2 * J J4+1 ) = 2 **D2*E 


00008160 

0594 


GO TO 102 


00008170 


c 



00008180 


c 

MATRIX TERMS AT AN INTERFACE FOR J=INF AND 

I BETWEEN FSW1 AND FSW2 

00008190 


c 



00008200 

0595 


222 XK « FLOAT ( K ) 


00008210 

0596 


A(KJltJJl) = 3.*<C55(M)+C55( P) ) 


00008220 

0597 


A(KJ1 t JJ6) = -4.*C55(P) 


00008230 

0598 


A(KJ1 t JJ8) = — 4.*C55 ( M ) 


00008240 

0599 


A ( K J 1 » J J 10 ) a C 5 5 ( M ) 


00008250 

0600 


A(KJ1,JJ12) a C551P) 


00008260 


c 



00008270 

0601 


AUJltJJl+1) = 3. 3 MC45(M)+C45(P) ) 


00008280 

0602 


A ( K J l f J J6+ 1 ) * -4.*C45(P) 


00008290 

0603 


A I K J 1 » J J 8+ 1 ) = —4 . *C45 ( M ) 


00008300 

0604 


A ( KJ1 * JJ10+1 ) = C45 ( M ) 


00008310 

0605 


A(KJ1, JJ12+1 ) = C45 ( P ) 


00008320 


c 



00008330 

0606 


A { K J 1 » J J2 + 2 ) = < C45 { P ) -C45 ( M ) )/XK 


00008340 

0607 


A ( K J 1 » J J4+2 ) = (C45(M)-C45(P) )/XK 


00008350 


c 



00008360 

0608 


A(KQltJJl) = 3.*(C45(M)+C45< P) > 


00008370 

0609 


A ( KQ 1 » J J 6 ) = -4.*C45(P) 


00008380 

0610 


A(KQ1,JJ8) = -4,*C45(M) 


00008390 

0611 


A ( KQ1 t J J10 ) = C45 ( M ) 


00008400 

0612 


A(KQ1, JJ12) = C45 ( P) 


00008410 


c 



00008420 

0613 


A(K01 f JJl+l) = 3. * ( C44 ( M ) +C44 ( P ) ) 


00008430 

0614 


A (KQ1 t J J6+1 ) = -4.*C44(P) 


00008440 

0615 


A ( KOI » J J8+ 1 ) = — 4. *C44 { M ) 


00008450 

0616 


A ( KQ1 , J J10+ 1 ) = C44 { M ) 


00008460 

0617 


A(KQ1, JJ12+1 J = C44 { P ) 


00008470 


c 



00008480 

0618 


A { KQ 1 » J J2+2 ) = ( C44 ( P ) -C44 ( M ) ) /XK 


00008490 

0619 


A ( KQ1 t J J4+2 ) = ( C44 { M ) — C44 ( P ) ) /XK 


00008500 


c 



00008510 

0620 


A(KQ2,JJ2) = (C36( P)-C36(M) )/XK 


00008520 

0621 


A ( KQ2 f J J4 ) = (C36(M)-C36(P))/XK 


00008530 


c 



00008540 

0622 


AIKQ2, JJ2+1) a (C23( P)-C23(M) ) /XK 


00008550 

0623 


A t KQ2 * J J4 + 1 ) a (C23{M)-C23(P))/XK 


00008560 


c 



00008570 

0624 


A ( KQ2 f J J 1 + 2 ) a 3.’MC33(M)+C33(P) ) 


00008580 

0625 


A { KQ2 ♦ J J6+2 ) = ~4.*C33(P) 


00008590 

0626 


A ( KQ2 * J J8+2 ) = -4.*C33(M) 


00008600 

0627 


A ( KQ2 1 J J 10+2 ) = C33 ( M ) 


00008610 

0628 


A { KQ2 » J J 12 + 2 ) = C33 ( P ) 


00008620 

0629 


X ( J01 ) = 0. 


00008630 


c 



00008640 

0630 


CZ1 = (C13(P)-C13(M) )*C3 + (C23( P)-C23(M) ) *BV + (C36 ( P ) -C36 { M ) ) *81100008650 

0631 


CZ2 = ( C 1 3 ( P ) — C 1 3 ( M ) )*C2+(C23<P)-C23(M) ) *DV + 2 .* ( C36 { P ) -C36( M ) )*C4 

00008660 


c 



00008670 

0632 


X(JJ1) = 0. 


00008680 

0633 


X ( JQ2 ) a 2.*H*<CZ1 + CZ2*Z) 


00008690 

0634 


GO TO 102 


00008700 
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C 

C 

FREE SURFACE MATRIX TERMS AT ANY INTERFACE WHERE I=LAW, J=INF OR AT 

00008710 

00008720 


C 

THE FREE SURFACE POINT I=LAW, J=L AT 

00008730 

0635 

C 

223 A ( K J 1 f J J 1 ) = 3.*C66(M) 

00008740 

00008750 

0636 


A { K J 1 » J J2 ) = -4.*C66(M) 

00008760 

0637 


A ( KJ 1 , J J 13 ) = C66(M) 

00008770 

0638 

c 

A ( K J ItJJl + l) = 3.*C26(M) 

00008780 

00008790 

0644 

c 

A(KQltJJl) = 3.*C26(M) 

00008860 

00008870 

0645 


A ( KQ1 « J J2 ) = — 4. *C26 ( M ) 

00008880 

0646 


A ( KQ1 v J J 13 ) = C26 ( M ) 

00008890 

0647 

c 

A ( K 0 1 1 J J 1+ 1 ) = 3 • *C22 ( M ) 

00008900 

00008910 

0648 


A ( KQ1 * J J2+1 ) = — 4. *C22 ( M ) 

00008920 

0649 


A (KQ1 f JJ13+1 ) - C22 ( M ) 

00008930 

0650 

c 

A(KQl f JJl+2) - 3.*C23(M) 

00008940 

00008950 

0651 


A { KQ 1 1 J J 8+2 ) * -4. *C23 ( M) 

00008960 

0652 


A ( KQ1 * J J 10+2 ) * C23 ( M ) 

00008970 

0653 

c 

A ( KQ2 » J J 1 ) = 3.*C45(M) 

00008980 

00008990 

0654 


A ( KQ2* J J8 ) * -4.*C45(M) 

00009000 

0655 


A ( KQ2 ♦ J J 10 ) = C45 ( M ) 

00009010 

0656 

c 

A { KQ2 t J J 1 + 1 ) = 3.*C44(M) 

00009020 

00009030 

0657 


A I KQ2 ? JJ8+ 1 ) = -4.*C44(M) 

00009040 

0658 


A(KQ2 t JJ10+1) = C44 { M ) 

00009050 

0659 

c 

A { KQ2 r J Jl + 2 ) = 3.*C44(M) 

00009060 

00009070 

0660 


A ( KQ2 » J J2+2 ) = —4. *C44 ( M ) 

00009080 

0661 


A ( KQ2 * J J 1 3+2 ) = C44 ( M ) 

00009090 

0662 

c 

CY1 = C12(M)*C3 + C22IM)*BV + C26<M)*BU 

00009100 

00009110 

■0663 


CY2 = C12(M)*C2 + C 22 ( M ) *DV + 2.*C26( M ) *C4 

00009120 

0664 


CXY 1 = C16(M)*C3 + C26(M)*BV + C66(M)*BU 

00009130 

0665 


CXY2 = C16(M)*C2 + C26(M)*DV + 2. *C66( M ) *C4 

00009140 

0666 

c 

X(JJl) = -2.*H*(CXY1 + CXY2*Z) 

00009150 

00009160 

0667 


X(JQ1) = -2.*H*(CY1 + C Y2*Z ) 

00009170 

0668 


X ( JQ2 ) = 0. 

00009180 

0669 


GO TO 102 

00009190 


c 

c 

MATRIX TERMS TO FIX THE RIGID TRANSLATIONS 

00009200 

00009210 

0670 

c 

203 A(KJlfJJl) = 1.0 

00009220 

00009230 

0671 


A ( KQ 1 » J J 1+ 1 ) = 1.0 

00009240 

0672 


A(KQ2 T JJl+2) = 1.0 

00009250 

0673 

c 

X( JJ1 ) = 0. 

00009260 

00009270 

0674 


X(JQ1) = 0. 

00009280 
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0675 

0676 

0677 

0678 


0679 

0680 
0681 

0682 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 


0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 


0711 

0712 


X ( JQ2 ) = 0. 

GO TO 102 
C 

202 IFC I.EQ.l) GO TO 220 
I F ( I.EQ.LAW) GO TO 223 
C 

C FREE SURFACE MATRIX TERMS FOR I BETWEEN 1 AND LAW AND J=LAT 
C 

A(KJl'JJl) = 3.*C55(M) 

A(KJlfJJ8) = -4.*C55(M> 

A ( K J 1 * J J 10 ) - C55 ( M ) 

C 

A(KJ1,JJ1+1) = 3.*C45(M) 

A(KJl t JJ8+l) = — 4.*C45(M) 

AtKJl, JJ10+1) = C45 ( M ) 

C 

A(KQltJJl) = 3 • *C45 { M ) 

A ( KQ1 1 J J 8 ) = -4.*C45(M) 

A(KQltJJlO) = C45 ( M ) 

C 

A(KQ1,JJ1+1) = 3.*C44(M) 

A(KQ1,JJ8+1) = — 4 »*C44 ( M ) 

A ( KOI » J J 10+ 1 ) = C44 ( M ) 

C 

A ( KQ2 t J J 1+2 ) = 3 . *C 33 ( M) 

A ( KQ2 ♦ J J 8+2 ) = -4.*C33(M) 

A ( KQ2 ♦ J J 10+2 ) = C33 { M ) 

C 

CZ1 = C13(M)*C3 + C23(M)*BV + C36(M)*BU 
CZ2 = C 1 3 ( M ) #C2 + C23(M)*DV + 2 . *C36 ( M ) *C4 
C 

X( JJ1 ) = 0. 

X( J01) = 0. 

XIJQ2) = -2.*H*(CZ1 + CZ2*Z) 

c 

I F ( I.E0.FSW1) GO TO 231 
I F ( I.E0.FSW2) GO TO 231 
IF ( I .GT.FSW1 .AND. I .LT.FSW2 ) GO TO 234 


C 

C IF 
C 


C 


C 


C 


I IS BETWEEN 1 AND FSW1 OR BETWEEN FSW2 AND LAW 


A ( KJ 1 r J J2+2 ) = — C45 ( M ) 
A(KJl,JJ4+2) = C45 ( M ) 


A(K01»JJ2+2) = -C44 ( M ) 
A(KQl»JJ4+2) = C44 ( M ) 


A(KQ2,JJ2) = -C36 ( M ) 
A { KQ2 * J J4 ) = C36( M) 


A ( K02 » J J2+ 1 ) = -C23 ( M) 
A ( K02 * J J4+1 ) = C 23 ( M ) 
GO TO 102 


C 

C CASE WHERE I=FSW1 OR FSW2 AND J= L A T 
C 

231 XK = FLOAT(K) 

D 1 = 2.*(XK-1.)/XK 


CONTINUE BELOW 


00009290 

00009300 

00009310 

00009320 

00009330 

00009340 

00009350 

00009360 

00009370 

00009380 

00009390 

00009400 

00009410 

00009420 

00009430 

00009440 

00009450 

00009460 

00009470 

00009480 

00009490 

00009500 

00009510 

00009520 

00009530 

00009540 

00009550 

00009560 

00009570 

00009580 

00009590 

00009600 

00009610 

00009620 

00009630 

00009640 

00009650 

00009660 

00009670 

00009680 

00009690 

00009700 

00009710 

00009720 

00009730 

00009740 

00009750 

00009760 

00009770 

00009780 

00009790 

00009800 

00009810 

00009820 

00009830 

00009840 

00009850 

00009860 
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0713 



D2 = 2.*XK/ { XK+1 . ) 

00009870 

0714 



D3 = 2./ ( < XK+1. ) *XK) 

00009880 


C 




00009890 

0715 



I F ( I.EQ.FSW2) GO TO 232 

00009900 


C 




00009910 

0716 



AtKJl, JJ1+2) 

» D 1*C 45 ( M ) 

00009920 

0717 



AIKJl , JJ2+2) 

* ~D2*C45(M) 

00009930 

0718 



A< KJ1, JJ4+2) 

* D3*C45 ( M ) 

00009940 


c 




00009950 

0719 



A ( KQ1 , JJ1+2 ) 

s D1*C44(M) 

00009960 

0720 



A(KQ1 t JJ2+2) 

» -*D2*C44(M) 

00009970 

0721 



A ( KQ1 , J J4+2 ) 

* D3*C44(M) 

00009980 


c 




00009990 

0722 



A ( KQ2 y J J 1 ) = 

D1*C36{ M) 

00010000 

0723 



A(KQ2,JJ2) = 

-D2*C36( M) 

00010010 

0724 



A ( KQ2 f J J4 ) = 

D3*C36< M) 

00010020 


c 




00010030 

0725 



A ( KQ2* JJ1+1 ) 

= D1*C23(M) 

00010040 

0726 



A(KQ2 , J J2+1 ) 

* -D2*C23(M) 

00010050 

0727 



A ( KQ2 * J J4+ 1 ) 

= D3*C23(M) 

00010060 

0728 



GO TO 102 


00010070 


c 




00010080 

0729 


232 

A ( K J l y J J 1 + 2 ) 

= -D1*C45(M) 

00010090 

0730 



A ( K J 1 y J J2+2 ) 

= -D3*C45(M> 

00010100 

0731 



A ( K J 1 y J J4+2 ) 

= D2*C45{M) 

00010110 


c 




00010120 

0732 



A ( KOI , J J 1+2 ) 

= -Dl*C44(M) 

00010130 

0733 



A( KOI, JJ2+2) 

* -D3*C44(M) 

00010140 

0734 



A ( KQ1 f J J4+2 ) 

= D2*C44(M) 

00010150 


c 




00010160 

0735 



A ( KQ2 y J J 1 ) = 

— D 1 *C36 ( M ) 

00010170 

0736 



A ( KQ2 , J J 2 ) - 

-D3*C36(M) 

00010180 

0737 



A ( KQ2 y J J4 ) = 

D2*C36(M ) 

00010190 


c 




00010200 

0738 



A ( KQ2 , J J 1+ 1 ) 

= -D1*C23(M) 

00010210 

0739 



AIKQ2, JJ2+1 ) 

= -D3*C23(M) 

00010220 

0740 



A(KQ2, JJ4+1 ) 

= D2*C23 { M ) 

00010230 

0741 



GO TO 102 


00010240 


c 




00010250 


c 

CASE 

: WHERE I IS 

BETWEEN FSW1 AND FSW2 AND J=L AT 

00010260 


c 




00010270 

0742 


234 

XK = FLOAT ( K 

) 

00010280 

0743 



A ( K J 1 y J J2 + 2 ) 

= -C45(M)/XK 

00010290 

0744 



AtKJly JJ4+2) 

= C45 { M ) /XK 

00010300 


c 




00010310 

0745 



AIKQl, JJ2+2) 

= -C44 { M ) / XK 

00010320 

0746 



AIKOl y JJ4+2) 

= C44(M)/XK 

00010330 


c 




00010340 

0747 



A ( KQ2 y J J2 ) = 

-C36(M) /XK 

00010350 

0748 



A ( KQ2 y J J4 ) = 

C 36 ( M ) /XK 

00010360 


c 




00010370 

0749 



A(KQ2y JJ2+1) 

= -C 23 ( M ) /XK 

00010380 

0750 



A ( KQ2 y J J4+ 1 } 

= C 2 3 ( M ) /XK 

00010390 


c 




00010400 

0751 


102 

CONTINUE 


00010410 


c 




00010420 


c 

FORM THE NONSYMETRIC BANDED MATRIX AX 

00010430 


c 




00010440 
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0752 

IL = K J 1+3* ( NODE-1 ) 


00010450 

0753 

IN = IL+2 
C 


00010460 

00010470 

0754 

DO 103 IK = IL , IN 


00010480 

0755 

II = IK-IL+1 
C 


00010490 

00010500 

0756 

DO 104 JK= 1 * NB AND 


00010510 

0757 

JJ = IK+JK-IBW-1 


00010520 

0758 

I F ( IK.LE.IBW1) JJ = JK 


00010530 

0759 

IF! JJ.GT. JQMAX) GO TO 105 


00010540 

0760 

AX! JK, IK) = A ( II , JJ) 


000105.50 

0761 

GO TO 104 


00010560 

0762 

105 AX( JK, IK) = 0.0 


00010570 

0763 

104 CONTINUE 


00010580 

0764 

103 CONTINUE 


00010590 

0765 

101 CONTINUE 


00010600 

0766 

100 CONTINUE 
C 


00010610 

00010620 

0767 

REWIND 9 


00010630 

0768 

WRITE (9) (( AX( J, I) » J= 1 , NB AND ) ,1=1, JQMAX) 


00010640 

0769 

WR I TE ( 9 ) (X( I) ,1-1, JQMAX) 


00010650 

0770 

END FILE 9 


00010660 

0771 

REWINO 9 
C 


00010670 

00010680 

0772 

NBD - NBAND+1 


00010690 

0773 

00 107 1=1, JQMAX 


00010700 

0774 

AX ( NBD , I ) = X! I ) 


00010710 

0775 

107 CONTINUE 
C 

C WRITE! 6,4000 ) 


00010720 

00010730 

00010740 


C4000 FORMAT ( 1H1 , ' EQUATION', 35X, 'THE BANDED MATRIX 
C CALL RITE(1, JQMAX, NBD, JQMAX, NBD, AX) 

C WRITE! 6,4003) 

TERMS AX ( I » J ) 1 

// J00010750 
00010760 
00010770 


C4003 FORMAT ! 1H1 , 45X, •*** THE LOAD VECTOR X ( I ) *** * 
C WRITE! 6,4004) (X(I), 1=1, JQMAX) 

C4004 FORMAT! 28! 2X, 10D12.3 / )) 

C 

III ) 

00010780 

00010790 

00010800 

00010810 

0776 

CALL TRMS TR ! AX , JQMAX, NBD , IBW, I0W, NBAND, 
C 

DT, RT, ET) 

00010820 

00010830 

0777 

WRITE ( 6 ,4006) ET, RT, DT 


00010840 

0778 

4006 FORMAT!/// • ERROR CONDITION OF SOLVER ROUTINE 
1 'RANK IS ', F6.1, 5X , 'DETERMINANT = G10.3) 

IS ', F4.1, 5X , 

00010850 

00010860 

0779 

I F ( ET. EQ. 1 • ) STOP 1 
C 


00010870 

00010880 

0780 

00 108 1=1, JQMAX 


00010890 

0781 

X!I) = AX! 1,1 ) 


00010900 

0782 

108 CONTINUE 
C 


00010910 

00010920 

0783 

READ! 9 ) ( ( AX ( J , I ) , J=l, NBAND) , 1=1, JQMAX) 


00010930 

0784 

RE AD ! 9 ) (R( I ) , 1=1, JQMAX) 

C 
C 


00010940 

00010950 

00010960 


C i***********#*************#**^********************* *******£**********00010970 
C 00010980 
C OUTPUT OF THE NODAL DISPLACEMENTS, U» V, W 00010990 
C 00011000 


C **$*:*:*>!£** a*************** *$$**$* ##**$** *****$>!£*** a******* c****#**>}c$#**:fcfc0001 10 10 

c 00011020 
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0785 


WRITE ( 6 t 650 ) 

00011030 

0786 


J = 1 

00011040 

0787 


DO 12 IK = 1, LAW 

00011050 

0788 


DO 11 JK = 1, LAT 

00011060 

0789 


WR I TE ( 6,651 ) J, X ! 3* J— 2 ) * X(3*J-1>, X!3*J) 

00011070 

0790 


J = J+l 

00011080 

0791 

11 

CONTINUE 

00011090 

0792 


WRITE! 6,653 ) 

00011100 

0793 

12 

C 

CONTINUE 

00011110 

00011120 

0794 


WR ITE ! 6 , 9950 ) 

00011130 

0795 

9950 

FORMAT! 1H1, 5X, 'EQUATION', 5X, '*** THE ACCURACY TEST, TEST-RtI) 
1 ***», 10X, •*** THE AVERAGE ABSOLUTE ERROR **** ///) 

00011140 

00011150 

0796 


ERR = 0. DO 

00011160 

0797 


DO 9990 1 = 1 , JQMAX 

00011170 

0798 


TEST = 0. DO 

00011180 

0799 


DO 9960 J = 1 »NB AND 

00011190 

0800 


IC = I+J-IBW-1 

00011200 

0801 


IFd.LE.IBWU IC - J 

00011210 

0802 


IF! IC.GT. JQMAX) GO TO 9970 

00011220 

0803 


TEST = TEST+AX! J, I)*X( IC) 

00011230 

0804 

9960 

CONTINUE 

00011240 

0805 

9970 

TEST = TEST-R(I) 

00011250 

0806 


ERR - ERR+DABS! TEST) 

00011260 

0807 


AV6 = ERR/I 

00011270 

0808 


WR I TE ( 6, 9980 ) I, TEST, AVE 

00011280 

0809 

9980 

FORMAT ! 5X , I4,10X, G15.8, 32X, G15.8) 

00011290 

0810 

9990 

C 

CONTINUE 

00011300 

00011310 


C ********** ***************************** *********** *********** ********* 0001 1320 
c 00011330 


c 

c 

CALCULATION OF THE STRAIN (S) AND STRESS (T) 

00011340 

00011350 


C ********************************************** ************************QQQ1^3£Q 
C 00011370 

0811 


SXM = SXMAX * 1.E06 

00011380 

0812 


SXE = C3E * 1.E06 

00011390 

0813 


WR ITE! 6,670) SXM, SXE 

00011400 

0814 


WRITE(6,671) 

00011410 

0815 


HR = l./(2.*H) 

00011420 

0816 

c 

XK = FLOAT(K) 

00011430 

00011440 

0817 


DO 399 1*1, LAW 

00011450 

0818 

c 

DO 398 J*l, LAT 

00011460 

00011470 

0819 


11=1-1 

00011480 

0820 


12=1-2 

00011490 

0821 


NODE = LAT*I1+J 

00011500 

0822 


JJ1 = 3*( LAT*Il+J)-2 

00011510 

0823 


JJ2 = 3*<LAT*I2+J)-2 

00011520 

0824 


JJ3 = 3* ( L AT* I 2+ J ) — 5 

00011530 

0825 


JJ4 = 3* ( L AT* 1+ J ) — 2 

0001 1540 

0826 


JJ5 = 3*!LAT*I+J)+1 

00011550 

0827 


JJ6 = 3* ! LAT* I 1+ J ) + 1 

00011560 

0828 


JJ7 = 3* ( LAT* I2+J ) +1 

00011570 

0829 


JJ8 = 3* ( LAT* I 1+ J ) -5 

0001 1580 

0830 


JJ9 = 3* t LAT* 1 + J ) —5 

00011590 

0831 


JJ10 =3*(LAT*Il+J)-8 

0001 1600 
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0832 



JJ11 = 3* ( LAT* ( 1+ 1 ) + J ) -2 


00011610 

0833 



JJ12 =3*UAT*Il + J)+4 


00011620 

0834 



JJ13 -3*(LAT*( I-3)+J)-2 


00011630 


C 




00011640 

0835 



Z = ( F LOAT t J ) - ( FLOAT ( L AT ) + 1 * ) / 2. ) *H 


00011650 

0836 



SX * C2*Z + C3 


00011660 


C 




00011670 

0837 



IF(I.EQ.l) GD TO 385 


00011680 

0838 



I F { I.EQ.LAW) GO TO 386 


00011690 

0839 



IF( I.GT.FSW1.AND.I.LT.FSW2) GO TO 382 


00011700 

0840 



IF( I,E0,FSW1.0R*I,EQ*FSW2) GO TO 383 


00011710 


C 




00011720 

0841 



HI = H 


00011730 

0842 



H2 = HI 


00011740 

0843 



GO TO 384 


00011750 


C 




00011760 

0844 


382 

HI = XK*H 


00011770 

0845 



H2 = HI 


00011780 

0846 



GO TO 384 


00011790 


C 




00011800 

0847 


383 

HI = H 


00011810 

0848 



H2 = XK*H 


00011820 

0849 



I F ( I.EQ.FSW1) GO TO 384 


00011830 

0850 



HI = XK*H 


00011840 

0851 



H2 = H 


00011850 


C 




00011860 

0852 


384 

H 1 2 = H1/H2 


00011870 

0853 



H21 = H2/H1 


00011880 

0854 



HRO = (H2-H1)/(H1*H2> 


00011890 

0855 



HRS = 1 • / ( H 1+H2 ) 


00011900 


C 




00011910 

0856 



SY = HRS* ( H12*X ( JJ4+1 ) -H2 1*X ( JJ2+1 ) ) +HRD*X { JJ1+1 ) + DV*Z 

+ BV 

00011920 

0857 



SX Y = HRS*<H12*X( JJ4)-H21*X( JJ2) ) + HRD*X(JJ1) + 2.*C4*Z 

+ BU 

0001 1930 

0858 



SYZ I = HRS*(H12*X( J J4+2 ) -H21*X( JJ2+2 ) )+HRO*X( JJ1+2 ) 


00011940 

0859 



GO TO 387 


0001 1950 


c 




00011960 

0860 


385 

SY = HR*(4.*X( JJ4+1)-3.*X{ JJ1+1)— X< JJ11+1) ) + DV*Z + BV 


00011970 

0861 



SXY = HR*(4.*X( JJ4)-3.*X( JJ1 ) -X ( JJ11 > ) + 2.*C4*Z + BU 


00011980 

0862 



SYZ I = HR*(4.*X( JJ4+2)-3,*X( JJl+2)-X( JJ11+2) ) 


00011990 

0863 



GO TO 387 


00012000 


c 




00012010 

0864 


386 

SY = HR*(3.*Xl JJ1+1J+X(JJ13+1)-4,*X(JJ2+1}) + DV*Z + BV 


00012020 

0865 



SXY = HR*(3.*X( JJ1)+X( JJ13)-4,*X( JJ2) ) + 2.*C4*Z + BU 


00012030 

0866 



SYZ I = HR*(3.*XUJl+2)+X{ JJ13+2)-4.*X( JJ2 + 2) ) 


00012040 


c 




00012050 

0867 


387 

DO 392 M=l t NL AY 


00012060 

0868 



IF(M.EQ. l.AND. J.GT. INF{ 1 ) ) GO TO 392 


00012070 

0869 



IFtM.EO.l) GO TO 388 


00012080 

0870 



I F { J.LE. INF (M-l) .OR, J,GT, INF(M) ) GO TO 392 


00012090 

0871 


388 

IF(J.EQ.l) GO TO 389 


00012100 

0872 



IF( J.EQ. INF{M) .OR.J.EQ.LAT) GO TO 390 


00012110 


c 




00012120 

0873 



SZ = HR*(X{ JJ6+2)-X{ JJ8+2) ) 


00012130 

0874 



SYZ J = HR*(X( JJ6+1)-X( JJ8+1) ) 


00012140 

0875 



SXZ = HR* < X( JJ6) — X( JJ8) ) 


00012150 

0876 



GO TO 391 


00012160 


c 




00012170 

0877 


389 

SZ = HR*(4.*X( JJ6+2)-3.*X( JJl+2)-X( JJ12+2)) 


00012180 
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0878 



SYZJ = HR*(4.*X( JJ6+1 >-3.*X( JJ1 + 1 )-X( JJ12+1 ) ) 

00012190 

0879 



SXZ - HR*(4.*X( JJ6)-3.*X( JJ1)-X< JJ12) ) 

00012200 

0880 



GO TO 391 

00012210 


C 



00012220 

0881 


390 

SZ = HR*( 3.*X( JJl+2)+X< JJ10+2)-4.*X< JJ8+2) > 

00012230 

0882 



SYZJ = HR*(3.*X( JJ1+1)+X( JJ10+1)-4.*X( JJ8+1) ) 

00012240 

0883 



SXZ = HR*(3.*X< JJ1)+X( JJ10)-4.*X( JJ8J ) 

00012250 


c 



00012260 

0884 


391 

SYZ - SYZI + SYZJ 

00012270 


c 



00012280 


c 


CALCULATION OF THE STRESS (T) 

00012290 


c 



00012300 

0885 



TX = C11(M)*SX + C12{M)*SY + C13(M)*SZ + C16(M)*SXY 

00012310 

0886 



TY = C 12< M) *SX + C22(M)*SY + C23(M)*SZ + C26(M)*SXY 

00012320 

0887 



TZ * C13 { M ) *SX + C23(M)*SY + C33(M)*SZ + C36(M)*SXY 

00012330 


c 



00012340 

0888 



TYZ = C44(M)*SYZ + C45(M)*SXZ 

00012350 

0889 



TXZ = C45 ( M ) *SYZ + C55(M)*SXZ 

00012360 

0890 



TXY * C16 ( M ) *SX + C26{M)*SY + C36<M)*SZ + C66(M)*SXY 

00012370 


c 



00012380 

0891 



WR ITE (6, 672 ) NODE, TX, TY, TZ, TYZ, TXZ, TXY, SY, SZ, SYZ, SXZ , SXYOOO 12390 

0892 



WR I TE ( 6t 397 ) SX 

00012400 


c 



00012410 


c 

STRESS AND STRAINS JUST ABOVE AN INTERFACE 

00012420 


c 



00012430 

0893 



IF( J.NE.INF(M).OR.J.EO.LAT) GO TO 392 

00012440 

0894 



P = M+l 

00012450 

0895 



SZ = HR*(4.*X( JJ6+2)-3.*X( JJl+2)-X( JJ12+2) ) 

00012460 

0896 



SYZJ = HR*(4.*X{ JJ6+ 1 >-3.*X( JJ1+1 )-X( JJ12+1 > ) 

00012470 

0897 



SXZ = HR*(4.*X( JJ6)-3.*X< JJ1 )-X( JJ12) ) 

00012480 

0898 



SYZ = SYZI + SYZJ 

00012490 


c 



00012500 

0899 



TX = C11(P)*SX + C12(P)*SY + C 1 3 ( P ) #S Z + C16(P)*SXY 

00012510 

0900 



TY = C12(P)*SX + C22 ( P ) *SY + C23(P)*SZ + C26(P)*SXY 

00012520 

0901 



TZ = C13(P)*SX + C23(P)*SY + C33(P)*SZ + C36(P)*SXY 

00012530 


c 



00012540 

0902 



TYZ = C44(P)*SYZ + C45(P)*SXZ 

00012550 

0903 



TXZ = C45(P)*SYZ + C55(P)*SXZ 

00012560 

0904 



TXY = C 1 6 ( P ) * SX + C26 ( P ) * SY + C36(P)*SZ + C66(P)*SXY 

00012570 


c 



00012580 

0905 



WRITE(6,672) NODE, TX, TY, TZ, TYZ, TXZ, TXY, SY, SZ, SYZ, SXZ , SXYOOO 12590 


c 



00012600 

0906 


392 

continue 

00012610 

0907 


398 

CONTINUE 

00012620 

0908 



WR I TE ( 6 » 652 ) 

00012630 

0909 


399 

CONTINUE 

00012640 

0910 

9000 

CONTINUE 

00012650 


c 



00012660 


c 



c 



00012680 


c 


FORMATS 

00012690 


c 



00012700 


c 

**** 


00012710 


c 



00012720 

0911 


397 

FORMAT { 14X, 1P1E11.3/) 

00012730 

0912 


600 

FORMAT ( 1 H 1 1 44X, 44H*** UNIFORM BENDING OF A LAMINATED PLATE ***) 

00012740 

0913 


601 

FORMAT ( 5110) 

00012750 


c 



00012760 
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0916 

0917 

0918 

0919 

0920 

0921 

0922 

0923 

0924 

0925 

0926 

0927 

0928 

0929 


0930 

0931 

0932 

0933 
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MAIN 


DATE = 75007 


08/16/07 


602 FORMAT ( //// 5X, 18H*** INPUT DATA *** /// 

1 18X, ‘NUMBER OF LAYERS IN CROSS SECTION* NLAY = ', 14 // 

2 18X, 'NUMBER OF NODES ON VERTICAL AXIS* LAT = ', 14 // 

3 18X, 'NUMBER OF NODES ON HORIZONTAL AXIS* LAW = 14 /// 

4 18X, 37HCHANGE IN MESH WIDTH (FSW1) AT I = , 14 // 

5 18X, 37HCHANGE IN MESH WIDTH * FSW2 ) AT I = , 14 // 

6 18X, 37HMESH WIDTH MAGNIFICATION FACTOR, K = , 14 /) 

603 F0RMAT(8G12.5) 

604 FORMAT ( 1H1, 55X, 21H*** MATERIAL DATA *** //// 2X, 5HLAYER , 7X, 

1 3HE11, 9X, 3HE22, 9X, 3HE33, 9X, 3HE12, 9X, 3HE13, 9X, 

2 3HE23, 8X* 4HNU12, 4X, 4HNU13, 4X, 4HNU23 // ) 

605 FORMAT ( 3X* 12, 6X, 2PE10.3, 2(2X, 1PE10.3), 3(2X, 0PE10.3), 

1 313X, F5.2) / ) 


606 FORMAT ( 10G10. 3 I 

607 FORMAT*/// 18X, 26HFINE SIMULATION WIDTH, H 

608 FORMAT*// 18X, 9HLAYER NO., 2X, 


, F8 • 5 ) 


5X , 17H INTERFACE AT J = ,13) 

611 FORMAT*// 45X , 41H*** COEFFICIENTS OF THERMAL EXPANSION ***, /// 
1 IX, 5HLAYER, 8X, 5HTHETA, 12X, 3HAL1, 12X, 3HAL2, 12X, 


3HAL3, 12X, 3HAL6, 12X, 4HAL1P, 
/// ) 


11X, 4HAL2P, 11X, 4HAL3P 


00012770 
00012780 
00012790 
00012800 
00012810 
00012820 
00012830 
00012840 
00012850 
00012860 
00012870 
00012880 
00012890 
00012900 
00012910 
00012920 
00012930 
00012940 
00012950 
00012960 
00012970 
00012980 
00012990 
00013000 
00013010 
00013020 
00013030 
00013040 
00013050 
00013060 
00013070 
00013080 
00013090 
00013100 
00013110 
00013120 
00013130 
00013140 
00013150 

650 FORMAT* 1H1 // 10X, • *** GRID POINT DISPLACEMENT FUNCTIONS ***' ///00013160 

1 16X, 5H NODE, 5X, 14HU-D I SPLACEMENT , 6X, 14HV-D I SPLAC EMENT , 0001 3 170 

2 6X , 14HW-DISPL ACEMENT /// ) 00013180 

00013190 

651 FORMATdOX, 110, 3E20.6 // ) 00013200 

652 FORMAT*// 12H ********** // ) 00013210 

653 FORMAT*// 10X, 12H ********** // ) 00013220 

00013230 

670 FORMAT* 1H1, 10X, 77H*** OUTPUT STRESSES AND STRAINS FOR A MAXIMUM 00013240 
1L0NGITUDINAL BENDING STRAIN OF , F6.0, 22H MICRO -INCHES/INCH AND /00013250 

2 48X, 40H AN APPLIED AXIAL EXTENS I ON AL STRAIN OF , F6.0, 00013260 

3 19H M ICRO- INCHES/ INCH. // 10X, 'NOTE: INTERFACE NODES ARE REPEA T0001 3270 
4ED WITH VALUES GIVEN BELOW AND ABOVE THE INTERFACE RESPECTIVELY.' 00013280 

00013290 
00013300 
00013310 
000-13320 
00013330 
00013340 
00013350 
00013360 
00013370 
00013380 
00013390 


613 FORMAT*// 53X, 26H*** STIFFNESS MATRICES *** /// IX, 

1 11HLAYER /THETA, 21X, 12HX-Y-Z MATRIX, 44X, 

2 18HX-Y-Z PRIME MATRIX /// ) 

614 FORMAT * 2X, 12, 9X, F5.1, 5X, 7*5X, E10.3)) 

620 FORMAT * 2X, 12, 5X, 1P12E10.3 // 19X, 5E10.3, 10X, 5E10.3 // 29X, 

1 4E10.3, 20X, 4E10.3 // 1X,0PF5.1, 33X, 1P3E10.3, 30X, 

2 3E10.3 // 49X, 2E10.3, 40X, 2E10.3 // 59X, E10.3, 50X, 

3 E10.3 /// ) 


5 //// ) 

671 FORMAT* IX, 5HN0DE , 5X, 5HSIG-X, 6X, 5HSIG-Y, 6X, 5HSIG-Z, 6X, 

1 6HTAU-YZ, 5X, 6HTAU-XZ , 5X, 6HTAU-XY, 5X, 5HEPS-Y, 6X, 

2 5HEPS-Z, 6X , 6HEPS-YZ, 5X, 6HEPS-XZ, 5X, 6HEPS-XY / 17X, 


3 5HEPS-X /// ) 

'672 FORMAT * IX, 13, 4X, 1P11E11.3 /) 

STOP 

END 
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0001 


RELEASE 2.0 


MATCON 


DATE = 75082 


19/49/20 


SUBROUTINE MATCON 


00013400 
00013410 

* * *i« »}< # ifc## #$$!<£ ifcjjojtsjtifcjic*# $$$$ $$### $4# $ $ $ $ 4; * * $ $ $ * $ *000 13420 

00013430 

CALCULATION OF LAMINATE LOAD CONSTANTS FOR A FULL CROSS SECTION 00013440 

00013450 

**********.*****.*#****** ************** ******** *************************000 13460 

00013470 

THIS SUBROUTINE IS GOOD FOR BENDING OF AN ARBITRARILY LAID UP 00013471 



C 

LAMINATE WHICH IS SYMMETRIC OR NONSYMMETRIC ABOUT THE MIDPLANE. 

00013472 


c 


00013480 


C 

THE CONSTANTS ARE C2 = INVERSE BENDING RADIUS 

00013490 


c 

C3E * APPLIED UNIFORM EXTENS ION AL STRAIN 

00013500 


c 

C3 = EXTENS I ONAL COUPLING DUE TO BENDING PLUS C3E 

00013510 


c 

C4 = IN-PLANE SHEAR COUPLING 

00013520 


c 


00013530 


c 

BU OCCURS IN THE FGTN. U(Y,Z) 

00013540 


c 

BV AND DV OCCUR IN THE FCTN. V(Y,Z) 

00013550 


c 


00013560 


c 

SXMAX (EFFECTIVELY THE LOAD INPUT) IS A MAXIMUM STRAIN 

00013570 


c 


00013580 

0002 


INTEGER ORDER 

00013590 


c 


00013600 

0003 


COMMON /MC/ C11(6),C12(6),C16(6),C22(6),C26(6),C66(6),C13(6), 

00013610 



1 C23 ( 6) ,C36 ( 6 ) ,C44(6) * C45 ( 6 > ,C5 5 ( 6 ) ,C33 ( 6 ) , AL 1 { 6 ) , AL2 { 6 ) , 

00013620 



2 AL3( 6) , AL6( 6) , C2 , C3 , C3E , C4, BU , DU, BV , DV , H, SXMAX , NLA Y , INF { 6) 

00013630 


c 


00013640 

0004 


DIMENSION A ( 3 « 3 ) , B(3,3), 0(3,3), QM(3,3) 

00013650 


c 


00013660 

0005 


DOUBLE PRECISION A, B, D 

00013670 


c 


00013680 

0006 


ORDER = 3 

00013690 


c 


00013700 

0007 


LAY = INF ( 1 ) — 1 

00013710 

0008 


HL = H*FLOAT(LAY) 

00013720 

0009 


HL2 = HL*#2/2. 

00013730 

0010 


HL 3 = HL**3/3. 

00013740 

Q01 1 


RN = FLOAT(NLAY) 

00013750 

0012 


RN2 = RN**2 

00013760 


c 


00013770 

0013 


DO 20 1=1,3 

00013780 

0014 


DO 20 J= 1 » 3 

00013790 

0015 


A ( I , J ) = 0 .DO 

00013800 

0016 


B( I, J) = 0. DO 

00013810 

0017 


D ( I , J ) = O.DO 

00013820 

0018 


20 CONTINUE 

00013830 


c 


00013840 

0019 


DO 30 1=1,3 

00013850 

0020 


DO 30 J = 1 , 3 

00013860 

0021 


DO 30 M= 1 T NL A Y 

00013870 

0022 


QM (1,1) = Cll (M)-C13(M)*C13(M)/C33(M) 

00013880 

0023 


QM (1,2) = C12(M)-C13(M)*C23(M)/C33(M) 

00013890 

0024 


OM (1,3) = C16(M)-C13(M)*C36(M)/C33<M) 

00013900 

0025 


QM ( 2 , 1 ) = QM ( 1 , 2 ) 

00013910 

0026 


QM (2,2) = C22 ( M ) — C23 (M)=*C23(M)/C33(M) 

00013920 

0027 


QM ( 2 , 3 ) = C 2 6 ( M ) — C 23{M)=i t C36(M)/C33(M) 

00013930 

0028 


QM (3,1) = QM( 1,3) 

00013940 

0029 


QM (3,2) = QM ( 2 , 3 ) 

00013950 
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0030 


QM<3,3) = C66{ M)-C36{M)*C36(M) /C33<M) 

00013960 


C 


00013970 


C 

NOTE THAT THE SUBSCRIPT 3 IN QM REPLACES A 6 IN STANDARD NOTATION. 

00013980 


C 

THE SAME IS TRUE BELOW IN A(I,J), B(I,J), D{I,J), ETC. 

00013990 


C 


00014000 

0031 


Ml = 2*M— 1 

00014010 

0032 


M2 = 3*M*(M-1)+1 

00014020 


c 


00014030 

0033 


A { I , J ) = A(I,J) + HL*QM(I,J) 

00014040 

0034 


B { 1 1 J ) = B(I r J) + HL2*QM{ I , J ) - ( M 1-RN ) 

00014050 

0035 


0(1, JJ = D { I , J ) + HL3*QM( I, J)*(M2-1.5*RN*M1+.75*RN2) 

00014060 

0036 


30. CONTINUE 

00014070 


c 


00014080 


c 

INVERT (A). STORE IN (A). 

00014090 


c 


00014100 

0037 


CALL MATIN4 (A, ORDER) 

00014110 


c 


00014120 


c 

MULTIPLY (A) INVERSE * (B). STORE IN A. 

00014130 


c 


00014140 

0038 


CALL MAMULT ( A, B, ORDER, A) 

00014150 


c 


00014160 


c 

MULTIPLY (B) * { A ) INVERSE * (B). STORE IN B. 

00014170 


c 


00014180 

0039 


CALL MAMULT ( B , A , ORDER, B ) 

00014190 


c 


00014200 

0040 


DO 40 1=1,3 

00014210 

0041 


DO 40 J= 1 , 3 

00014220 

0042 


At I, J) = -l.*A( I , J ) 

00014230 

0043 


D ( I,J) = D ( I, J) - B( I, J) 

00014240 

0044 


40 CONTINUE 

00014250 


c 


00014260 


c 

INVERT NEW MATRIX (D). THE RESULT IS D-PRIME. STORE IN D. 

00014270 


c 


00014280 

0045 


CALL MAT I N4 (D, ORDER) 

00014290 


c 


00014300 


c 

MULTIPLY -(A) INVERSE * B * D-PRIME WHICH YIELDS B-PRIME. STORE IN B. 

00014310 


c 


00014320 

0046 


CALL MAMULT ( A , D , ORDE R , B ) 

00014330 


c 


00014340 


c 

DETERMINE THE LOAD CONSTANTS. MINUS C2 IMPLIES A SMILING PLATE. 

00014350 


c 


00014360 

0047 


ZMAX = RN*HL/2. 

00014370 

0048 


C2 = -D(1,1)*SXMAX/(B( 1,1) +D ( 1,1 ) *ZMAX ) 

00014380 

0049 


RATIO = C2/DI 1 , 1 ) 

00014390 


c 


00014400 

0050 


C3 = B< 1,1) *RAT 10 + C3E 

00014410 

0051 


C4 = . 5*D( 1,3)* RATIO 

00014420 

0052 


BU = B ( 3 , 1 ) *R AT 1 0 

00014430 

0053 


BV = B ( 2 , 1 ) *RAT IQ 

00014440 

0054 


DV = D(1,2)*RATI0 

00014450 


c 

RATIO = -RATIO 

00014460 

0055 


WRITE (6,50) 

00014470 

0056 


50 FORMAT (///// 48X, 35H*** THE LAMINATE LOAD CONSTANTS *** /// ) 

00014480 

0057 


WR I TE ( 6 , 60 ) C2, C3, C4, BU, BV, DV, RATIO 

00014490 

0058 


60 FORMAT { 1 C2 = *, 1PE10.3, 4X, *C3 = E10.3, 4X, 'C4 = ',E10.3, 

00014500 



1 4X, 1 BU = », E10.3, 4X , *BV = ', E10.3, 4X, * DV = *, E10.3, 

,00014510 



2 4X t «MT = * , E10.3 ) 

00014520 

0059 


RETURN 

00014530 

0060 


END 

00014540 
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MATCON 


DATE 


75082 


19/49/20 


FORTRAN IV 

Gl 

♦OPTIONS 

IN 

♦OPTIONS 

IN 

♦statist 

ICS* 

♦STATIST 

ICS* 


RELEASE 2.0 


EFFECT* NOTERM, NOIO, EBCDIC, SOURCE, NOLI ST, NODECK, LOAD, NOMAP, NOTEST 
EFFECT* NAME = MATCON , L INECNT = 60 

SOURCE STATEMENTS = 60, PROGRAM SIZE = 2060 

NO DIAGNOSTICS GENERATED 


FORTRAN IV Gl RELEASE 2.0 


MAMULT 


DATE = 75082 


19/49/20 


OOOl 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 
0013 


C 

C 

C 

C 


SUBROUTINE MAMULT < B , C , N , A ) 

MAMULT POSTMULTIPLIES MATRIX (B) BY MATRIX (C) AND STORES THE 
RESULT IN MATRIX (A) WHERE N IS THE ORDER OF THE MATRICES. 


DOUBLE PRECISION A,B,C,$UM 
DIMENSION A ( N ,N ) , B(N,N), C(N,N) 
DO 1 1=1, N 
DO 1 J= 1 , N 
SUM = 0. 

DO 2 K= 1 , N 

SUM = SUM + B t I,KJ*C(K,J) 

2 CONTINUE 
At I , J) = SUM 

1 continue 

RETURN 

END 


00014550 

00014551 

00014552 

00014553 

00014554 

00014560 

00014570 

00014580 

00014590 

00014600 

00014610 

00014620 

00014630 

00014640 

00014650 

00014660 

00014670 


FORTRAN IV Gl RELEASE 2.0 


MAMULT 


DATE = 75082 


19/49/20 


*OPT IONS IN 
♦OPTIONS IN 
♦STATISTICS* 
♦STATISTICS* 


EFFECT* NOTERM, NO ID, EBCDIC, SOURC E , NOL I ST , NODECK , LOAD 
EFFECT* NAME = MAMULT , L INECNT = 60 

SOURCE STATEMENTS = 13, PROGRAM SIZE = 

NO DIAGNOSTICS GENERATED 


, NOMAP, NOTES 
702 


T 


FORTRAN IV Gl RELEASE 2.0 MAT IN4 DATE = 75082 19/49/20 


OOOl 



subroutine matiN4 ( array, n) 


00014680 


C 




00014681 


C 

MATIN4 INVERTS THE MATRIX (ARRAY) WHICH ! 

IS OF ORDER N. 

00014682 


C 




00014683 

0002 



DIMENSION ARRAY ( N , N ) 


00014690 

0003 



DOUBLE PRECISION ARRAY 


00014700 

0004 



DO 604 1=1, N 


00014710 

0005 



STORE = ARRAY ( 1,1) 


00014720 

0006 



ARRAY t I, I) = 1. 


00014730 

0007 



DO 601 J=l,N 


00014740 

0008 


601 

ARRAY ( I , J ) = ARRAYt I , J) /STORE 


00014750 

0009 



DO 604 K=1,N 


00014760 

0010 



I F ( K — 1 ) 602,604,602 


00014770 

0011 


602 

STORE = ARRAY ( K , I ) 


00014780 

0012 



ARRAY ( K , I ) *= 0. 


00014790 

0013 



DO 603 J= 1 , N 


00014800 

0014 


603 

ARRAY ( K , J ) = ARRA Y ( K , J ) - STORE*ARRAY( 

I , J ) 

00014810 

0015 


604 

CONTINUE 


00014820 

0016 



RETURN 


00014830 

0017 



END 


00014840 
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TRMSTR 


DATE 


75007 


08/16/07 


0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 


SUBROUTINE TRMSTR ( A ,N, ND , NL D, NRD, NED, D, R, E > 

C 

C TRMSTR IS THE SUBROUTINE TRIMSS WITH MATRIX A TRANSPOSED. 

C THE SIMULTANEOUS SOLUTIONS IS GAUSSIAN ELIMINATION, 

C MODIFIED TO TAKE ADVANTAGE OF THE REDUCED MATRIX. THE 

C ROUTINE ALSO USES PARTIAL PIVOTING TO REDUCE ROUNDOFF ERROR. 

C INPUT 

C 1 A 

C 
C 
C 
C 
C 
C 
C 
c 
c 
c 
c 
c 
c 
c 
c 

c 2 N 

C 3 ND 

C 

C 4 NLD 

C 
C 

C 5 NRD 

C 
C 

C 6 NED 

C OUTPUT 

C 1 A 

C 

C 2D 

C 3 R 

C 4 E 

C 
C 
C 
C 

C SUBROUTINE TRMSTR ( A, N ,ND, NLD , NRD , NED, D , R, E ) 

DIMENSION A ( ND » 1 ) 

DOUBLE PRECISION A,D,Y,W,S 
XI = 1. 

LI = 1 
E = 0. 

R = 0. 

D= 1 . 

ND1=NED+1 
M = NLD 
NM1=N-1 
DO 1 1=1, NM1 
I F ( I.GT.(N-NLD) )M=M-1 
NN= I+M-l 
DO 2 1 1 = I , NN 


FIRST LOCATION OF COEFFICIENT MATRIX, I. E. A(l,l>. 
THE BAND ELEMENTS IN EACH ROW MUST BE LEFT 
JUSTIFIED AND EXTEND TO THE RIGHT M PLACES 
(M=MIN( N,NLD+NRD+1 ) • IF IN ANY PARTICULAR ROW 
THERE ARE ONLY K BAND ELEMENTS AND K IS LESS 
THAN M, THEN THE M-K RIGHT MOST ELEMENTS OF THAT 
ROW WILL BE SET TO ZERO. THE ROW WHOSE LEFT 
MOST COLUMN IN THE FULL BLOWN MATRIX CONTAINS 
A NON-ZERO ELEMENT MUST BE THE FIRST ROW OF THE 
REDUCED MATRIX AND ETC. THE COLUMN TO THE 
IMMEDIATE RIGHT OF THE' REDUCED MATRIX (FORMED AS 
ABOVE) MUST CONTAIN THE RIGHT HAND SIDE OF THE 
EQUATION SET IN QUESTION. IT SHOULD NOW BE 
OBVIOUS THAT AN N X N+l FULL BLOWN SYSTEM WOULD 
BE REDUCED BY THE ABOVE METHOD TO AN N X M+l 
system. 

NUMBER OF SIMULTANEOUS EQUATIONS TO BE SOLVED. 
VARIABLE DIMENSION INTEGER. MUST BE EQUAL TO 
ROW DIMENSION OF A IN CALLING PROGRAM. 

MAXIMUM NUMBER OF BAND ELEMENTS TO THE LEFT 
OF PRINCIPAL DIAGONAL IN ANY ROW OF SYSTEM TO 
BE DETERMINED. 

MAXIMUM NUMBER OF BAND ELEMENTS TO THE RIGHT 
OF PRINCIPAL DIAGONAL IN ANY ROW OF SYSTEM TO 
BE DETERMINED. 

NED=MIN ( N, NLD+NRD+1 ) 

THE FIRST COLUMN OF A CONTAINS THE SOLUTION 
VECTOR. 

CONTAINS DETERMINANT OF A. 

CONTAINS RANK OF A. 

E=0. t SOLUTION O.K. E=l., A SINGULAR. 

E=2 . » SOLUTION ATTEMPTED, BUT A ILL CONDITIONED 
OR SINGULAR. IN THIS CASE SOLUTIONS SHOULD BE 
CHECKED TO ASSURE VALIDITY. 


00014850 

00014860 

00014870 

00014880 

00014890 

00014900 

00014910 

00014920 

00014930 

00014940 

00014950 

00014960 

00014970 

00014980 

00014990 

00015000 

00015010 

00015020 

00015030 

00015040 

00015050 

00015060 

00015070 

00015080 

00015090 

00015100 

00015110 

00015120 

00015130 

00015140 

00015150 

00015160 

00015170 

00015180 

00015190 

00015200 

00015210 

00015220 

00015230 

00015240 

00015250 

00015260 

00015270 

00015280 

00015290 

00015300 

00015310 

00015320 

00015330 

00015340 

00015350 

00015360 

00015370 

00015380 

00015390 

00015400 

00015410 

00015420 
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0016 



IF(DABS< A( 1, I 1 ) .GE.DABSfAl 1,1 1 + 1 ) ) ) GO TO 2 

00015430 

0017 



0 = -D 


00015440 

0018 



DO 3 1 , ND1 


00015450 

0019 



Y= A ( J , I) 


0001546*0 

0020 



A(J, I )=A( J,II+1) 


00015470 

0021 


3 

A ( J 1 1 1 + 1 ) =Y 


00015480 

0022 


2 

CONTINUE 


00015490 



C 

D-D*A( 1, I) 


00015500 

0023 



IF ( A ( If 1 ) .EQ. 0.) GO TO 10 


00015510 

0024 



GO TO (5, 13) «L1 


00015520 

0025 


13 

IF t DABS ( DABS ( ( X 1-A ( 1 , I ) ) /XI) -1 . ) 

.LT.l.E-07) E=2. 

00015530 

0026 



XI = Ad, I ) 


00015540 

0027 


5 

R = R + 1. 


00015550 

002 8 



LI = 2 


00015560 

0029 



DO 4 J=2 » ND1 


00015570 

0030 


4 

A { J , I ) =A ( J, I )/ A(l, I) 


00015580 

0031 



K= 1 + 1 


00015590 

0032 



NN= I+M 


00015600 

0033 



DO 1 II=K,NN 


00015610 

0034 



W-A{1, II ) 


00015620 

0035 



DO 6 J s 1 ,NED 


00015630 

0036 


6 

A( J, 1 1 )=A( J+l, I I ) —A < J+l, I )*W 


00015640 

0037 



A ( ND1 , I I ) = A ( NED , I I ) 


00015650 

0038 


1 

A< NED, I I ) =0 . 


00015660 

0039 



I F ( A ( 1 , N ) - EQ . 0 - ) GO TO 10 


00015670 

0040 



I F ( DABS ( DABS ( ( Xl-A( 1,N))/X1)-1.) 

• LT. 1 • E— 07 ) E = 2 . 

00015680 

0041 


9 

R = R + 1. 


00015690 

0042 



A ( 1 » N ) - A < ND 1 , N ) / A ( 1,N) 


00015700 

0043 



K-NM1 


00015710 

0044 



NN-2 


00015720 

0045 


8 

IF ( NN -GT -NED ) NN=NED 


00015730 

0046 



J = K+1 


00015740 

0047 



S=0. 


00015750 

0048 



DO 7 I =2 , NN 


00015760 

0049 



S=5+A ( 1 , J ) * A ( IfK) 


00015770 

0050 


7 

J = J+1 


00015780 

0051 



A(1 ,K) = A(ND1,K)-S 


00015790 

0052 



NN=NN+ 1 


00015800 

0053 



K=K-1 


00015810 

0054 



IF(K-NE-O) GO TO 8 


00015820 

0055 



RETURN 


00015830 

0056 


10 

E= 1 . 


00015840 

0057 



RETURN 


00015850 

0058 



END 


00015860 

FORTRAN 

IV 

G1 RELEASE 

2.0 TRMSTR 

DATE = 75007 08/16/07 


^OPTIONS 

IN EFFECT* 

NOTERM,NOID,EBCDIC,SOURCE,NOL 1ST 

, NODECK, LOAD, NOMAP, NOTE ST 


*OPT IONS 

IN EFFECT* 

NAME = TRMSTR , LINECNT = 

60 


*STATI 

ST 

ICS* SOURCE STATEMENTS = 58, PROGRAM 

SIZE = 2294 


*STATI 

STICS* NO DIAGNOSTICS GENERATED 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 


SUBROUTINE RITE! I DUM , NR , NC , MR, MC , A ) 
DOUBLE PRECISION A 
DIMENSION A(MR,MC) 

I PR INT= 12 

I F ( IDUM.NE. 1 ) I PR INT= 30 
I PR= IPRINT-1 
DO 35 K= 1 , NC , I PR I NT 
MAX= K+IPR 
IF(MAX.GT.NC) MAX=NC 
IF(K.NE.l) WR I TE ( 6 , 103 ) 

45 WRITE(6,102) (I,I=K,MAX) 

DO 40 J-l , NR 

40 WR I TE ( 6, 105 ) J , < A ( J , I ) , I=K , MAX ) 

35 CONTINUE 
RETURN 

101 FORMAT ( 6X , 3014) 

102 F0RMAT(6X, 12110) 

103 FORMAT( '1* ) 

104 FORMAT ( • ■ , 15,3014) 

105 FORMAT ( * *,15,12010.3) 

END 


00015870 

00015880 

00015890 

00015900 

00015910 

00015920 

00015930 

00015940 

00015950 

00015960 

00015970 

00015980 

00015990 

00016000 

00016010 

00016020 

00016030 

00016040 

00016050 

00016060 

00016070 


FORTRAN IV G1 RELEASE 2.0 


RITE 


DATE = 75007 08/16/07 


^OPTIONS IN EFFECT* N0T6RM, NO ID , EBCD IC , SOURCE ,NOL I ST , NODECK , LOAD, NOMAP , NOTEST 
^OPTIONS IN EFFECT* NAME = RITE , LINECNT = 60 

*STATISTICS* SOURCE STATEMENTS = 21, PROGRAM SIZE = 864 

*STAT I STICS* NO DIAGNOSTICS GENERATED 


*STATIST ICS* 


NO DIAGNOSTICS THIS STEP 
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